R Markdown

library(ggplot2)
library(tidyverse)
library(dplyr)
library(meta)
library(PRISMAstatement)
library(skimr)
library(MASS)
library(ggpubr)
library(patchwork)
### Site Level
final_data <- read.csv("Final Data.csv")
final_data_2022 <- final_data %>%
  filter(year == "2022")

shapiro.test(final_data_2022$animals)
## 
##  Shapiro-Wilk normality test
## 
## data:  final_data_2022$animals
## W = 0.76017, p-value = 7.037e-05
ggqqplot(final_data_2022$animals)

n1 <- glm(animals ~ shrub_density * aridity, family = "gaussian", data = final_data_2022)
n1
## 
## Call:  glm(formula = animals ~ shrub_density * aridity, family = "gaussian", 
##     data = final_data_2022)
## 
## Coefficients:
##           (Intercept)          shrub_density                aridity  
##              -113.435                 -6.222                 50.977  
## shrub_density:aridity  
##                 2.798  
## 
## Degrees of Freedom: 23 Total (i.e. Null);  20 Residual
## Null Deviance:       619900 
## Residual Deviance: 517700    AIC: 317.6
anova(n1, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: animals
## 
## Terms added sequentially (first to last)
## 
## 
##                       Df Deviance Resid. Df Resid. Dev Pr(>Chi)  
## NULL                                     23     619927           
## shrub_density          1    26978        22     592949  0.30732  
## aridity                1    71496        21     521453  0.09653 .
## shrub_density:aridity  1     3720        20     517733  0.70462  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(n1), "\n")
## Linear Model AIC: 317.6089
quadratic_model_1 <- glm(animals ~ shrub_density + I(shrub_density^2) + aridity + I(aridity^2) + shrub_density:aridity, family = "gaussian", data = final_data_2022)
quadratic_model_1
## 
## Call:  glm(formula = animals ~ shrub_density + I(shrub_density^2) + 
##     aridity + I(aridity^2) + shrub_density:aridity, family = "gaussian", 
##     data = final_data_2022)
## 
## Coefficients:
##           (Intercept)          shrub_density     I(shrub_density^2)  
##             -4704.343                 69.124                 -7.773  
##               aridity           I(aridity^2)  shrub_density:aridity  
##              2395.612               -287.024                  6.017  
## 
## Degrees of Freedom: 23 Total (i.e. Null);  18 Residual
## Null Deviance:       619900 
## Residual Deviance: 180300    AIC: 296.3
anova(quadratic_model_1, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: animals
## 
## Terms added sequentially (first to last)
## 
## 
##                       Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                     23     619927              
## shrub_density          1    26978        22     592949  0.100791    
## I(shrub_density^2)     1     2942        21     590007  0.587895    
## aridity                1    78131        20     511876  0.005227 ** 
## I(aridity^2)           1   315016        19     196860 2.052e-08 ***
## shrub_density:aridity  1    16537        18     180323  0.198864    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(quadratic_model_1), "\n")
## Quadratic Model AIC: 296.2959
if (AIC(quadratic_model_1) < AIC(n1)) {
  cat("Quadratic model is better.\n")
} else {
  cat("Linear model is better or equally good.\n")
}
## Quadratic model is better.
abund_dens_2022 <- ggplot(final_data_2022, aes(shrub_density, animals)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") +  labs(x = "Shrub Density per 20m Radius", y = "Animal Abundance") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
  theme(axis.title.x = element_blank()) + labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))

### Richness
n2 <- glm(richness ~ shrub_density * aridity, family = "gaussian", data = final_data_2022)
n2
## 
## Call:  glm(formula = richness ~ shrub_density * aridity, family = "gaussian", 
##     data = final_data_2022)
## 
## Coefficients:
##           (Intercept)          shrub_density                aridity  
##              -7.07856                0.25745                3.07545  
## shrub_density:aridity  
##              -0.04175  
## 
## Degrees of Freedom: 23 Total (i.e. Null);  20 Residual
## Null Deviance:       475 
## Residual Deviance: 334.1     AIC: 141.3
anova(n2, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: richness
## 
## Terms added sequentially (first to last)
## 
## 
##                       Df Deviance Resid. Df Resid. Dev Pr(>Chi)   
## NULL                                     23     474.96            
## shrub_density          1    6.431        22     468.53 0.534951   
## aridity                1  133.613        21     334.91 0.004681 **
## shrub_density:aridity  1    0.828        20     334.09 0.823795   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(n2), "\n")
## Linear Model AIC: 141.3094
quadratic_model_2 <- glm(richness ~ shrub_density + I(shrub_density^2) + aridity + I(aridity^2) + shrub_density:aridity, family = "gaussian", data = final_data_2022)
quadratic_model_2
## 
## Call:  glm(formula = richness ~ shrub_density + I(shrub_density^2) + 
##     aridity + I(aridity^2) + shrub_density:aridity, family = "gaussian", 
##     data = final_data_2022)
## 
## Coefficients:
##           (Intercept)          shrub_density     I(shrub_density^2)  
##            -128.52243               -0.17675                0.03871  
##               aridity           I(aridity^2)  shrub_density:aridity  
##              65.21718               -7.60162               -0.05383  
## 
## Degrees of Freedom: 23 Total (i.e. Null);  18 Residual
## Null Deviance:       475 
## Residual Deviance: 26.87     AIC: 84.82
anova(quadratic_model_2, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: richness
## 
## Terms added sequentially (first to last)
## 
## 
##                       Df Deviance Resid. Df Resid. Dev Pr(>Chi)    
## NULL                                     23     474.96             
## shrub_density          1    6.431        22     468.53  0.03794 *  
## I(shrub_density^2)     1  106.430        21     362.10  < 2e-16 ***
## aridity                1  101.019        20     261.08  < 2e-16 ***
## I(aridity^2)           1  232.882        19      28.20  < 2e-16 ***
## shrub_density:aridity  1    1.324        18      26.87  0.34636    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(quadratic_model_2), "\n")
## Quadratic Model AIC: 84.82182
if (AIC(quadratic_model_2) < AIC(n2)) {
  cat("Quadratic model is better.\n")
} else {
  cat("Linear model is better or equally good.\n")
}
## Quadratic model is better.
rich_dens_2022 <- ggplot(final_data_2022, aes(shrub_density, richness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") +  labs(x = "Shrub Density per 20m Radius", y = "Richness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
  theme(axis.title.x = element_blank()) + labs(tag = "B") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))

### Evenness
n3 <- glm(Average_Evenness ~ shrub_density * aridity, data = final_data_2022)
n3
## 
## Call:  glm(formula = Average_Evenness ~ shrub_density * aridity, data = final_data_2022)
## 
## Coefficients:
##           (Intercept)          shrub_density                aridity  
##             -0.059619              -0.002519               0.035851  
## shrub_density:aridity  
##              0.001281  
## 
## Degrees of Freedom: 23 Total (i.e. Null);  20 Residual
## Null Deviance:       0.1455 
## Residual Deviance: 0.1074    AIC: -51.7
anova(n3, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: Average_Evenness
## 
## Terms added sequentially (first to last)
## 
## 
##                       Df  Deviance Resid. Df Resid. Dev Pr(>Chi)  
## NULL                                      23    0.14552           
## shrub_density          1 0.0073326        22    0.13819  0.24267  
## aridity                1 0.0299710        21    0.10822  0.01817 *
## shrub_density:aridity  1 0.0007791        20    0.10744  0.70332  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(n3), "\n")
## Linear Model AIC: -51.70484
quadratic_model_3 <- glm(Average_Evenness ~ shrub_density + I(shrub_density^2) + aridity + I(aridity^2) + shrub_density:aridity, family = "gaussian", data = final_data_2022)
quadratic_model_3
## 
## Call:  glm(formula = Average_Evenness ~ shrub_density + I(shrub_density^2) + 
##     aridity + I(aridity^2) + shrub_density:aridity, family = "gaussian", 
##     data = final_data_2022)
## 
## Coefficients:
##           (Intercept)          shrub_density     I(shrub_density^2)  
##            -1.8496731             -0.0253566              0.0022251  
##               aridity           I(aridity^2)  shrub_density:aridity  
##             0.9526101             -0.1121067              0.0004442  
## 
## Degrees of Freedom: 23 Total (i.e. Null);  18 Residual
## Null Deviance:       0.1455 
## Residual Deviance: 0.01326   AIC: -97.91
anova(quadratic_model_3, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: Average_Evenness
## 
## Terms added sequentially (first to last)
## 
## 
##                       Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                     23   0.145519              
## shrub_density          1 0.007333        22   0.138186  0.001608 ** 
## I(shrub_density^2)     1 0.056383        21   0.081804 < 2.2e-16 ***
## aridity                1 0.018828        20   0.062975 4.308e-07 ***
## I(aridity^2)           1 0.049621        19   0.013354 2.285e-16 ***
## shrub_density:aridity  1 0.000090        18   0.013264  0.726529    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(quadratic_model_2), "\n")
## Quadratic Model AIC: 84.82182
if (AIC(quadratic_model_2) < AIC(n2)) {
  cat("Quadratic model is better.\n")
} else {
  cat("Linear model is better or equally good.\n")
}
## Quadratic model is better.
even_dens_2022 <- ggplot(final_data_2022, aes(shrub_density, Average_Evenness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") +  labs(x = "Shrub Density per 20m Radius", y = "Evenness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+
   labs(tag = "C") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
even_dens_2022

plot <- abund_dens_2022/rich_dens_2022/even_dens_2022
plot

### Site Level
final_data_2023 <- final_data %>%
  filter(year == "2023")

x1 <- glm(animals ~ shrub_density * aridity, family = "gaussian", data = final_data_2023)
x1
## 
## Call:  glm(formula = animals ~ shrub_density * aridity, family = "gaussian", 
##     data = final_data_2023)
## 
## Coefficients:
##           (Intercept)          shrub_density                aridity  
##               3.96923               -0.08162                0.93546  
## shrub_density:aridity  
##               0.20460  
## 
## Degrees of Freedom: 21 Total (i.e. Null);  18 Residual
## Null Deviance:       12650 
## Residual Deviance: 8381  AIC: 203.2
anova(x1, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: animals
## 
## Terms added sequentially (first to last)
## 
## 
##                       Df Deviance Resid. Df Resid. Dev Pr(>Chi)  
## NULL                                     21    12648.8           
## shrub_density          1   2257.7        20    10391.1  0.02766 *
## aridity                1   1570.4        19     8820.6  0.06628 .
## shrub_density:aridity  1    439.8        18     8380.8  0.33110  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(x1), "\n")
## Linear Model AIC: 203.1718
quadratic_model_x1 <- glm(animals ~ shrub_density + I(shrub_density^2) + aridity + I(aridity^2) + shrub_density:aridity, family = "gaussian", data = final_data_2023)
quadratic_model_x1
## 
## Call:  glm(formula = animals ~ shrub_density + I(shrub_density^2) + 
##     aridity + I(aridity^2) + shrub_density:aridity, family = "gaussian", 
##     data = final_data_2023)
## 
## Coefficients:
##           (Intercept)          shrub_density     I(shrub_density^2)  
##               68.7465                -5.2639                 0.3959  
##               aridity           I(aridity^2)  shrub_density:aridity  
##              -17.6832                 1.0897                 0.2494  
## 
## Degrees of Freedom: 21 Total (i.e. Null);  16 Residual
## Null Deviance:       12650 
## Residual Deviance: 7045  AIC: 203.4
anova(quadratic_model_x1, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: animals
## 
## Terms added sequentially (first to last)
## 
## 
##                       Df Deviance Resid. Df Resid. Dev Pr(>Chi)  
## NULL                                     21    12648.8           
## shrub_density          1  2257.72        20    10391.1  0.02355 *
## I(shrub_density^2)     1    55.46        19    10335.6  0.72267  
## aridity                1  1622.60        18     8713.0  0.05490 .
## I(aridity^2)           1  1189.72        17     7523.3  0.10022  
## shrub_density:aridity  1   478.47        16     7044.8  0.29720  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(quadratic_model_x1), "\n")
## Quadratic Model AIC: 203.3514
if (AIC(quadratic_model_x1) < AIC(x1)) {
  cat("Quadratic model is better.\n")
} else {
  cat("Linear model is better or equally good.\n")
}
## Linear model is better or equally good.
abund_dens_2023 <- ggplot(final_data_2023, aes(shrub_density, animals)) +
geom_smooth(method = "lm", se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") +  labs(x = "Shrub Density per 20m Radius", y = "Animal Abundance") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
  theme(axis.title.x = element_blank()) + labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
abund_dens_2023

x2 <- glm(richness ~ shrub_density * aridity, data = final_data_2023)
x2
## 
## Call:  glm(formula = richness ~ shrub_density * aridity, data = final_data_2023)
## 
## Coefficients:
##           (Intercept)          shrub_density                aridity  
##               1.79659                0.26232                0.28474  
## shrub_density:aridity  
##              -0.00865  
## 
## Degrees of Freedom: 21 Total (i.e. Null);  18 Residual
## Null Deviance:       95.32 
## Residual Deviance: 52.71     AIC: 91.66
anova(x2, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: richness
## 
## Terms added sequentially (first to last)
## 
## 
##                       Df Deviance Resid. Df Resid. Dev Pr(>Chi)   
## NULL                                     21     95.318            
## shrub_density          1  25.4984        20     69.820 0.003169 **
## aridity                1  16.3249        19     53.495 0.018219 * 
## shrub_density:aridity  1   0.7861        18     52.709 0.604382   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(x2), "\n")
## Linear Model AIC: 91.65558
quadratic_model_x2 <- glm(richness ~ shrub_density + I(shrub_density^2) + aridity + I(aridity^2) + shrub_density:aridity, data = final_data_2023)
quadratic_model_x2
## 
## Call:  glm(formula = richness ~ shrub_density + I(shrub_density^2) + 
##     aridity + I(aridity^2) + shrub_density:aridity, data = final_data_2023)
## 
## Coefficients:
##           (Intercept)          shrub_density     I(shrub_density^2)  
##             2.7345590             -0.0007208              0.0196188  
##               aridity           I(aridity^2)  shrub_density:aridity  
##             0.0227789              0.0152190             -0.0050956  
## 
## Degrees of Freedom: 21 Total (i.e. Null);  16 Residual
## Null Deviance:       95.32 
## Residual Deviance: 51.78     AIC: 95.26
anova(quadratic_model_x2, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: richness
## 
## Terms added sequentially (first to last)
## 
## 
##                       Df Deviance Resid. Df Resid. Dev Pr(>Chi)   
## NULL                                     21     95.318            
## shrub_density          1  25.4984        20     69.820  0.00500 **
## I(shrub_density^2)     1   0.4031        19     69.417  0.72413   
## aridity                1  17.3499        18     52.067  0.02059 * 
## I(aridity^2)           1   0.0894        17     51.977  0.86800   
## shrub_density:aridity  1   0.1998        16     51.778  0.80379   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(quadratic_model_x2), "\n")
## Quadratic Model AIC: 95.26343
if (AIC(quadratic_model_x2) < AIC(x2)) {
  cat("Quadratic model is better.\n")
} else {
  cat("Linear model is better or equally good.\n")
}
## Linear model is better or equally good.
rich_dens_2023 <- ggplot(final_data_2023, aes(shrub_density, richness)) +
geom_smooth(method = "lm", se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") +  labs(x = "Shrub Density per 20m Radius", y = "Richness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
  theme(axis.title.x = element_blank()) + labs(tag = "B") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
rich_dens_2023

x3 <- glm(Average_Evenness ~ shrub_density * aridity, data = final_data_2023)
x3
## 
## Call:  glm(formula = Average_Evenness ~ shrub_density * aridity, data = final_data_2023)
## 
## Coefficients:
##           (Intercept)          shrub_density                aridity  
##            -0.0315882              0.0055497              0.0122207  
## shrub_density:aridity  
##            -0.0001543  
## 
## Degrees of Freedom: 21 Total (i.e. Null);  18 Residual
## Null Deviance:       0.1608 
## Residual Deviance: 0.1076    AIC: -44.6
anova(x3, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: Average_Evenness
## 
## Terms added sequentially (first to last)
## 
## 
##                       Df Deviance Resid. Df Resid. Dev Pr(>Chi)  
## NULL                                     21    0.16080           
## shrub_density          1 0.013715        20    0.14708  0.12993  
## aridity                1 0.039182        19    0.10790  0.01048 *
## shrub_density:aridity  1 0.000250        18    0.10765  0.83796  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(x3), "\n")
## Linear Model AIC: -44.60481
quadratic_model_x3 <- glm(Average_Evenness ~ shrub_density + I(shrub_density^2) + aridity + I(aridity^2) + shrub_density:aridity, data = final_data_2023)
quadratic_model_x3
## 
## Call:  glm(formula = Average_Evenness ~ shrub_density + I(shrub_density^2) + 
##     aridity + I(aridity^2) + shrub_density:aridity, data = final_data_2023)
## 
## Coefficients:
##           (Intercept)          shrub_density     I(shrub_density^2)  
##            -7.236e-01              3.997e-03             -2.517e-05  
##               aridity           I(aridity^2)  shrub_density:aridity  
##             2.134e-01             -1.181e-02              2.469e-04  
## 
## Degrees of Freedom: 21 Total (i.e. Null);  16 Residual
## Null Deviance:       0.1608 
## Residual Deviance: 0.006624  AIC: -101.9
anova(quadratic_model_x3, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: Average_Evenness
## 
## Terms added sequentially (first to last)
## 
## 
##                       Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                     21   0.160797              
## shrub_density          1 0.013715        20   0.147082 8.624e-09 ***
## I(shrub_density^2)     1 0.020159        19   0.126922 2.991e-12 ***
## aridity                1 0.023110        18   0.103812 7.927e-14 ***
## I(aridity^2)           1 0.096719        17   0.007093 < 2.2e-16 ***
## shrub_density:aridity  1 0.000469        16   0.006624    0.2872    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(quadratic_model_x3), "\n")
## Quadratic Model AIC: -101.9453
if (AIC(quadratic_model_x1) < AIC(x1)) {
  cat("Quadratic model is better.\n")
} else {
  cat("Linear model is better or equally good.\n")
}
## Linear model is better or equally good.
even_dens_2023 <- ggplot(final_data_2023, aes(shrub_density, Average_Evenness)) +
geom_smooth(method = "lm", se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") +  labs(x = "Shrub Density per 20m Radius", y = "Evenness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+
   labs(tag = "C") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
even_dens_2023

plot2 <- abund_dens_2023/rich_dens_2023/even_dens_2023
plot2

### Temperature 2022 (Microsite/plot level)
t1 <- glm(animals ~ shrub_density * mean_temp, family = "gaussian", data = final_data_2022)
t1
## 
## Call:  glm(formula = animals ~ shrub_density * mean_temp, family = "gaussian", 
##     data = final_data_2022)
## 
## Coefficients:
##             (Intercept)            shrub_density                mean_temp  
##                641.3669                  22.3664                 -20.8447  
## shrub_density:mean_temp  
##                 -0.6221  
## 
## Degrees of Freedom: 23 Total (i.e. Null);  20 Residual
## Null Deviance:       619900 
## Residual Deviance: 389100    AIC: 310.8
anova(t1, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: animals
## 
## Terms added sequentially (first to last)
## 
## 
##                         Df Deviance Resid. Df Resid. Dev Pr(>Chi)   
## NULL                                       23     619927            
## shrub_density            1    26978        22     592949 0.238985   
## mean_temp                1   199846        21     393103 0.001351 **
## shrub_density:mean_temp  1     3970        20     389134 0.651496   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(t1), "\n")
## Linear Model AIC: 310.756
quadratic_temp_1 <- glm(animals ~ shrub_density + I(shrub_density^2) + mean_temp + I(mean_temp^2) + shrub_density:mean_temp, family = "gaussian", data = final_data_2022)
quadratic_temp_1
## 
## Call:  glm(formula = animals ~ shrub_density + I(shrub_density^2) + 
##     mean_temp + I(mean_temp^2) + shrub_density:mean_temp, family = "gaussian", 
##     data = final_data_2022)
## 
## Coefficients:
##             (Intercept)            shrub_density       I(shrub_density^2)  
##                4363.574                  155.190                   -6.801  
##               mean_temp           I(mean_temp^2)  shrub_density:mean_temp  
##                -305.916                    5.317                   -2.618  
## 
## Degrees of Freedom: 23 Total (i.e. Null);  18 Residual
## Null Deviance:       619900 
## Residual Deviance: 301200    AIC: 308.6
anova(quadratic_temp_1, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: animals
## 
## Terms added sequentially (first to last)
## 
## 
##                         Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                       23     619927              
## shrub_density            1    26978        22     592949 0.2042037    
## I(shrub_density^2)       1     2942        21     590007 0.6750263    
## mean_temp                1   240180        20     349827 0.0001516 ***
## I(mean_temp^2)           1     2772        19     347055 0.6840117    
## shrub_density:mean_temp  1    45820        18     301235 0.0979909 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(quadratic_temp_1), "\n")
## Quadratic Model AIC: 308.6113
if (AIC(quadratic_temp_1) < AIC(t1)) {
  cat("Quadratic model is better.\n")
} else {
  cat("Linear model is better or equally good.\n")
}
## Quadratic model is better.
abund_temp_2022 <- ggplot(final_data_2022, aes(mean_temp, animals)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") +  labs(x = "Mean Temperature (°C)", y = "Animal Abundance") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
  theme(axis.title.x = element_blank()) + labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
abund_temp_2022

t2 <- glm(richness ~ shrub_density * mean_temp, family = "gaussian", data = final_data_2022)
t2
## 
## Call:  glm(formula = richness ~ shrub_density * mean_temp, family = "gaussian", 
##     data = final_data_2022)
## 
## Coefficients:
##             (Intercept)            shrub_density                mean_temp  
##               32.305966                -0.039811                -1.022221  
## shrub_density:mean_temp  
##                0.005855  
## 
## Degrees of Freedom: 23 Total (i.e. Null);  20 Residual
## Null Deviance:       475 
## Residual Deviance: 120.1     AIC: 116.8
anova(t2, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: richness
## 
## Terms added sequentially (first to last)
## 
## 
##                         Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                       23     474.96              
## shrub_density            1     6.43        22     468.53    0.3008    
## mean_temp                1   348.04        21     120.49 2.703e-14 ***
## shrub_density:mean_temp  1     0.35        20     120.14    0.8088    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(t2), "\n")
## Linear Model AIC: 116.7632
quadratic_temp_2 <- glm(richness ~ shrub_density + I(shrub_density^2) + mean_temp + I(mean_temp^2) + shrub_density:mean_temp, family = "gaussian", data = final_data_2022)
quadratic_temp_2
## 
## Call:  glm(formula = richness ~ shrub_density + I(shrub_density^2) + 
##     mean_temp + I(mean_temp^2) + shrub_density:mean_temp, family = "gaussian", 
##     data = final_data_2022)
## 
## Coefficients:
##             (Intercept)            shrub_density       I(shrub_density^2)  
##              139.969166                -1.240892                 0.104613  
##               mean_temp           I(mean_temp^2)  shrub_density:mean_temp  
##               -9.215104                 0.152255                 0.007732  
## 
## Degrees of Freedom: 23 Total (i.e. Null);  18 Residual
## Null Deviance:       475 
## Residual Deviance: 74.51     AIC: 109.3
anova(quadratic_temp_2, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: richness
## 
## Terms added sequentially (first to last)
## 
## 
##                         Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                       23     474.96              
## shrub_density            1    6.431        22     468.53   0.21261    
## I(shrub_density^2)       1  106.430        21     362.10 3.965e-07 ***
## mean_temp                1  262.720        20      99.38 1.631e-15 ***
## I(mean_temp^2)           1   24.467        19      74.91   0.01505 *  
## shrub_density:mean_temp  1    0.400        18      74.51   0.75602    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(quadratic_temp_2), "\n")
## Quadratic Model AIC: 109.2983
if (AIC(quadratic_temp_2) < AIC(t2)) {
  cat("Quadratic model is better.\n")
} else {
  cat("Linear model is better or equally good.\n")
}
## Quadratic model is better.
rich_temp_2022 <- ggplot(final_data_2022, aes(mean_temp, richness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") +  labs(x = "Mean Temperature (°C)", y = "Richness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
  theme(axis.title.x = element_blank()) + labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
rich_temp_2022

t3 <- glm(Average_Evenness ~ shrub_density * mean_temp, family = "gaussian", data = final_data_2022)
t3
## 
## Call:  glm(formula = Average_Evenness ~ shrub_density * mean_temp, family = "gaussian", 
##     data = final_data_2022)
## 
## Coefficients:
##             (Intercept)            shrub_density                mean_temp  
##               0.4144890                0.0186680               -0.0124555  
## shrub_density:mean_temp  
##              -0.0005946  
## 
## Degrees of Freedom: 23 Total (i.e. Null);  20 Residual
## Null Deviance:       0.1455 
## Residual Deviance: 0.05211   AIC: -69.07
anova(t3, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: Average_Evenness
## 
## Terms added sequentially (first to last)
## 
## 
##                         Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                       23   0.145519              
## shrub_density            1 0.007333        22   0.138186   0.09341 .  
## mean_temp                1 0.082454        21   0.055732 1.847e-08 ***
## shrub_density:mean_temp  1 0.003627        20   0.052105   0.23807    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(t3), "\n")
## Linear Model AIC: -69.07199
quadratic_temp_3 <- glm(Average_Evenness ~ shrub_density + I(shrub_density^2) + mean_temp + I(mean_temp^2) + shrub_density:mean_temp, family = "gaussian", data = final_data_2022)
quadratic_temp_3
## 
## Call:  glm(formula = Average_Evenness ~ shrub_density + I(shrub_density^2) + 
##     mean_temp + I(mean_temp^2) + shrub_density:mean_temp, family = "gaussian", 
##     data = final_data_2022)
## 
## Coefficients:
##             (Intercept)            shrub_density       I(shrub_density^2)  
##               1.8333330               -0.0250452                0.0030448  
##               mean_temp           I(mean_temp^2)  shrub_density:mean_temp  
##              -0.1201345                0.0019981               -0.0002403  
## 
## Degrees of Freedom: 23 Total (i.e. Null);  18 Residual
## Null Deviance:       0.1455 
## Residual Deviance: 0.02914   AIC: -79.02
anova(quadratic_temp_3, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: Average_Evenness
## 
## Terms added sequentially (first to last)
## 
## 
##                         Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                       23   0.145519              
## shrub_density            1 0.007333        22   0.138186   0.03332 *  
## I(shrub_density^2)       1 0.056383        21   0.081804 3.601e-09 ***
## mean_temp                1 0.049707        20   0.032096 3.004e-08 ***
## I(mean_temp^2)           1 0.002571        19   0.029526   0.20760    
## shrub_density:mean_temp  1 0.000386        18   0.029140   0.62536    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(quadratic_temp_3), "\n")
## Quadratic Model AIC: -79.02002
if (AIC(quadratic_temp_3) < AIC(t3)) {
  cat("Quadratic model is better.\n")
} else {
  cat("Linear model is better or equally good.\n")
}
## Quadratic model is better.
even_temp_2022 <- ggplot(final_data_2022, aes(mean_temp, Average_Evenness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") +  labs(x = "Mean Temperature (°C)", y = "Evenness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
  theme(axis.title.x = element_blank()) + labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
even_temp_2022

### Temperature 2023 (Microsite/plot level)
z1 <- glm(animals ~ shrub_density * mean_temp, family = "gaussian", data = final_data_2023)
z1
## 
## Call:  glm(formula = animals ~ shrub_density * mean_temp, family = "gaussian", 
##     data = final_data_2023)
## 
## Coefficients:
##             (Intercept)            shrub_density                mean_temp  
##                 18.0908                  12.3265                  -0.1743  
## shrub_density:mean_temp  
##                 -0.4532  
## 
## Degrees of Freedom: 21 Total (i.e. Null);  18 Residual
## Null Deviance:       12650 
## Residual Deviance: 8513  AIC: 203.5
anova(z1, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: animals
## 
## Terms added sequentially (first to last)
## 
## 
##                         Df Deviance Resid. Df Resid. Dev Pr(>Chi)  
## NULL                                       21    12648.8           
## shrub_density            1  2257.72        20    10391.1  0.02890 *
## mean_temp                1  1313.77        19     9077.3  0.09558 .
## shrub_density:mean_temp  1   563.83        18     8513.5  0.27491  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(z1), "\n")
## Linear Model AIC: 203.5172
quadratic_temp_z1 <- glm(animals ~ shrub_density + I(shrub_density^2) + mean_temp + I(mean_temp^2) + shrub_density:mean_temp, family = "gaussian", data = final_data_2023)
quadratic_temp_z1
## 
## Call:  glm(formula = animals ~ shrub_density + I(shrub_density^2) + 
##     mean_temp + I(mean_temp^2) + shrub_density:mean_temp, family = "gaussian", 
##     data = final_data_2023)
## 
## Coefficients:
##             (Intercept)            shrub_density       I(shrub_density^2)  
##              1703.91419                  4.63774                 -0.02144  
##               mean_temp           I(mean_temp^2)  shrub_density:mean_temp  
##              -139.96488                  2.88050                 -0.14516  
## 
## Degrees of Freedom: 21 Total (i.e. Null);  16 Residual
## Null Deviance:       12650 
## Residual Deviance: 5706  AIC: 198.7
anova(quadratic_temp_z1, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: animals
## 
## Terms added sequentially (first to last)
## 
## 
##                         Df Deviance Resid. Df Resid. Dev Pr(>Chi)   
## NULL                                       21    12648.8            
## shrub_density            1  2257.72        20    10391.1 0.011869 * 
## I(shrub_density^2)       1    55.46        19    10335.6 0.693345   
## mean_temp                1  1517.38        18     8818.2 0.039146 * 
## I(mean_temp^2)           1  3064.70        17     5753.5 0.003375 **
## shrub_density:mean_temp  1    47.11        16     5706.4 0.716289   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(quadratic_temp_z1), "\n")
## Quadratic Model AIC: 198.716
if (AIC(quadratic_temp_z1) < AIC(z1)) {
  cat("Quadratic model is better.\n")
} else {
  cat("Linear model is better or equally good.\n")
}
## Quadratic model is better.
abund_temp_2023 <- ggplot(final_data_2023, aes(mean_temp, animals)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") +  labs(x = "Mean Temperature (°C)", y = "Animal Abundance") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
  theme(axis.title.x = element_blank()) + labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
abund_temp_2023

z2 <- glm(richness ~ shrub_density * mean_temp, family = "gaussian", data = final_data_2023)
z2
## 
## Call:  glm(formula = richness ~ shrub_density * mean_temp, family = "gaussian", 
##     data = final_data_2023)
## 
## Coefficients:
##             (Intercept)            shrub_density                mean_temp  
##                -3.72064                  0.97821                  0.33694  
## shrub_density:mean_temp  
##                -0.03281  
## 
## Degrees of Freedom: 21 Total (i.e. Null);  18 Residual
## Null Deviance:       95.32 
## Residual Deviance: 66    AIC: 96.6
anova(z2, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: richness
## 
## Terms added sequentially (first to last)
## 
## 
##                         Df Deviance Resid. Df Resid. Dev Pr(>Chi)   
## NULL                                       21     95.318            
## shrub_density            1  25.4984        20     69.820 0.008364 **
## mean_temp                1   0.8608        19     68.959 0.628021   
## shrub_density:mean_temp  1   2.9560        18     66.003 0.369258   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(z2), "\n")
## Linear Model AIC: 96.60373
quadratic_temp_z2 <- glm(richness ~ shrub_density + I(shrub_density^2) + mean_temp + I(mean_temp^2) + shrub_density:mean_temp, family = "gaussian", data = final_data_2023)
quadratic_temp_z2
## 
## Call:  glm(formula = richness ~ shrub_density + I(shrub_density^2) + 
##     mean_temp + I(mean_temp^2) + shrub_density:mean_temp, family = "gaussian", 
##     data = final_data_2023)
## 
## Coefficients:
##             (Intercept)            shrub_density       I(shrub_density^2)  
##              66.1097226                0.6326525                0.0007925  
##               mean_temp           I(mean_temp^2)  shrub_density:mean_temp  
##              -5.4527083                0.1192890               -0.0197340  
## 
## Degrees of Freedom: 21 Total (i.e. Null);  16 Residual
## Null Deviance:       95.32 
## Residual Deviance: 61.32     AIC: 98.99
anova(quadratic_temp_z2, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: richness
## 
## Terms added sequentially (first to last)
## 
## 
##                         Df Deviance Resid. Df Resid. Dev Pr(>Chi)   
## NULL                                       21     95.318            
## shrub_density            1  25.4984        20     69.820   0.0099 **
## I(shrub_density^2)       1   0.4031        19     69.417   0.7457   
## mean_temp                1   0.6475        18     68.769   0.6811   
## I(mean_temp^2)           1   6.5745        17     62.195   0.1903   
## shrub_density:mean_temp  1   0.8706        16     61.324   0.6337   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(quadratic_temp_z2), "\n")
## Quadratic Model AIC: 98.98617
if (AIC(quadratic_temp_z2) < AIC(z2)) {
  cat("Quadratic model is better.\n")
} else {
  cat("Linear model is better or equally good.\n")
}
## Linear model is better or equally good.
rich_temp_2023 <- ggplot(final_data_2023, aes(mean_temp, richness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") +  labs(x = "Mean Temperature (°C)", y = "Richness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
  theme(axis.title.x = element_blank()) + labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
rich_temp_2023

z3 <- glm(Average_Evenness ~ shrub_density * mean_temp, family = "gaussian", data = final_data_2023)
z3
## 
## Call:  glm(formula = Average_Evenness ~ shrub_density * mean_temp, family = "gaussian", 
##     data = final_data_2023)
## 
## Coefficients:
##             (Intercept)            shrub_density                mean_temp  
##               -0.879560                 0.033490                 0.039643  
## shrub_density:mean_temp  
##               -0.001127  
## 
## Degrees of Freedom: 21 Total (i.e. Null);  18 Residual
## Null Deviance:       0.1608 
## Residual Deviance: 0.04301   AIC: -64.79
anova(z3, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: Average_Evenness
## 
## Terms added sequentially (first to last)
## 
## 
##                         Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                       21   0.160797              
## shrub_density            1 0.013715        20   0.147082   0.01659 *  
## mean_temp                1 0.100583        19   0.046499 8.713e-11 ***
## shrub_density:mean_temp  1 0.003485        18   0.043014   0.22718    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(z3), "\n")
## Linear Model AIC: -64.78675
quadratic_temp_z3 <- glm(Average_Evenness ~ shrub_density + I(shrub_density^2) + mean_temp + I(mean_temp^2) + shrub_density:mean_temp, family = "gaussian", data = final_data_2023)
quadratic_temp_z3
## 
## Call:  glm(formula = Average_Evenness ~ shrub_density + I(shrub_density^2) + 
##     mean_temp + I(mean_temp^2) + shrub_density:mean_temp, family = "gaussian", 
##     data = final_data_2023)
## 
## Coefficients:
##             (Intercept)            shrub_density       I(shrub_density^2)  
##               5.642e+00                4.705e-03               -1.426e-04  
##               mean_temp           I(mean_temp^2)  shrub_density:mean_temp  
##              -5.012e-01                1.114e-02                5.355e-05  
## 
## Degrees of Freedom: 21 Total (i.e. Null);  16 Residual
## Null Deviance:       0.1608 
## Residual Deviance: 0.0005175     AIC: -158
anova(quadratic_temp_z3, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: Average_Evenness
## 
## Terms added sequentially (first to last)
## 
## 
##                         Df Deviance Resid. Df Resid. Dev Pr(>Chi)    
## NULL                                       21   0.160797             
## shrub_density            1 0.013715        20   0.147082   <2e-16 ***
## I(shrub_density^2)       1 0.020159        19   0.126922   <2e-16 ***
## mean_temp                1 0.085588        18   0.041335   <2e-16 ***
## I(mean_temp^2)           1 0.040811        17   0.000524   <2e-16 ***
## shrub_density:mean_temp  1 0.000006        16   0.000518   0.6562    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(quadratic_temp_z3), "\n")
## Quadratic Model AIC: -158.0316
if (AIC(quadratic_temp_z3) < AIC(z3)) {
  cat("Quadratic model is better.\n")
} else {
  cat("Linear model is better or equally good.\n")
}
## Quadratic model is better.
even_temp_2023 <- ggplot(final_data_2023, aes(mean_temp, Average_Evenness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") +  labs(x = "Mean Temperature (°C)", y = "Evenness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
  theme(axis.title.x = element_blank()) + labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
even_temp_2023

### Test colinearity 2022
your_model <- glm(animals ~ shrub_density + I(shrub_density^2) + mean_temp + I(mean_temp^2) + aridity + I(aridity^2) + mean_humidity + I(mean_humidity^2), family = poisson(link = "log"), data = final_data_2022)
your_model
## 
## Call:  glm(formula = animals ~ shrub_density + I(shrub_density^2) + 
##     mean_temp + I(mean_temp^2) + aridity + I(aridity^2) + mean_humidity + 
##     I(mean_humidity^2), family = poisson(link = "log"), data = final_data_2022)
## 
## Coefficients:
##        (Intercept)       shrub_density  I(shrub_density^2)           mean_temp  
##           68.06110             0.63570            -0.05244             6.70852  
##     I(mean_temp^2)             aridity        I(aridity^2)       mean_humidity  
##           -0.14576           -74.57567             7.80189             2.59059  
## I(mean_humidity^2)  
##           -0.04495  
## 
## Degrees of Freedom: 23 Total (i.e. Null);  15 Residual
## Null Deviance:       4036 
## Residual Deviance: 416.3     AIC: 562.1
summary(your_model)
## 
## Call:
## glm(formula = animals ~ shrub_density + I(shrub_density^2) + 
##     mean_temp + I(mean_temp^2) + aridity + I(aridity^2) + mean_humidity + 
##     I(mean_humidity^2), family = poisson(link = "log"), data = final_data_2022)
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         68.06110   30.34372   2.243 0.024897 *  
## shrub_density        0.63570    0.02693  23.603  < 2e-16 ***
## I(shrub_density^2)  -0.05244    0.00238 -22.037  < 2e-16 ***
## mean_temp            6.70852    1.96996   3.405 0.000661 ***
## I(mean_temp^2)      -0.14576    0.04251  -3.429 0.000607 ***
## aridity            -74.57567   21.89055  -3.407 0.000657 ***
## I(aridity^2)         7.80189    2.34090   3.333 0.000860 ***
## mean_humidity        2.59059    0.75099   3.450 0.000561 ***
## I(mean_humidity^2)  -0.04495    0.01281  -3.509 0.000450 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 4036.16  on 23  degrees of freedom
## Residual deviance:  416.33  on 15  degrees of freedom
## AIC: 562.14
## 
## Number of Fisher Scoring iterations: 6
your_model
## 
## Call:  glm(formula = animals ~ shrub_density + I(shrub_density^2) + 
##     mean_temp + I(mean_temp^2) + aridity + I(aridity^2) + mean_humidity + 
##     I(mean_humidity^2), family = poisson(link = "log"), data = final_data_2022)
## 
## Coefficients:
##        (Intercept)       shrub_density  I(shrub_density^2)           mean_temp  
##           68.06110             0.63570            -0.05244             6.70852  
##     I(mean_temp^2)             aridity        I(aridity^2)       mean_humidity  
##           -0.14576           -74.57567             7.80189             2.59059  
## I(mean_humidity^2)  
##           -0.04495  
## 
## Degrees of Freedom: 23 Total (i.e. Null);  15 Residual
## Null Deviance:       4036 
## Residual Deviance: 416.3     AIC: 562.1
a <- anova(your_model, test = "Chisq")
a
## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: animals
## 
## Terms added sequentially (first to last)
## 
## 
##                    Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                  23     4036.2              
## shrub_density       1   205.26        22     3830.9 < 2.2e-16 ***
## I(shrub_density^2)  1    22.68        21     3808.2  1.91e-06 ***
## mean_temp           1  2700.17        20     1108.0 < 2.2e-16 ***
## I(mean_temp^2)      1     0.00        19     1108.0 0.9503221    
## aridity             1   669.01        18      439.0 < 2.2e-16 ***
## I(aridity^2)        1     3.95        17      435.1 0.0468959 *  
## mean_humidity       1     6.48        16      428.6 0.0109107 *  
## I(mean_humidity^2)  1    12.27        15      416.3 0.0004594 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(car)
vif_results <- car::vif(your_model)
print(vif_results)
##      shrub_density I(shrub_density^2)          mean_temp     I(mean_temp^2) 
##           68.80459           72.66716        24431.18644        29383.10054 
##            aridity       I(aridity^2)      mean_humidity I(mean_humidity^2) 
##       235991.70543       224724.98700        15974.68265        10607.89891
cor_matrix <- cor(final_data_2022[, c("shrub_density", "mean_temp", "aridity", "mean_humidity")])
print(cor_matrix)
##               shrub_density   mean_temp     aridity mean_humidity
## shrub_density    1.00000000  0.02938781  0.02585627    0.04761062
## mean_temp        0.02938781  1.00000000 -0.85229523   -0.98132444
## aridity          0.02585627 -0.85229523  1.00000000    0.80272028
## mean_humidity    0.04761062 -0.98132444  0.80272028    1.00000000
library(lmerTest)
library(performance)
model <- lmer(animals ~ shrub_density + mean_temp + mean_humidity + (1|year), data = final_data)
check_collinearity(model)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##           Term  VIF    VIF 95% CI Increased SE Tolerance Tolerance 95% CI
##  shrub_density 1.11 [1.01,  2.77]         1.05      0.90     [0.36, 0.99]
## 
## Moderate Correlation
## 
##           Term  VIF    VIF 95% CI Increased SE Tolerance Tolerance 95% CI
##      mean_temp 7.50 [4.80, 12.13]         2.74      0.13     [0.08, 0.21]
##  mean_humidity 7.68 [4.90, 12.42]         2.77      0.13     [0.08, 0.20]
anova(model)
## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## shrub_density  26141   26141     1 41.354  2.4335 0.12639  
## mean_temp      61743   61743     1 41.634  5.7477 0.02108 *
## mean_humidity  12699   12699     1 39.750  1.1822 0.28347  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ranova(model)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## animals ~ shrub_density + mean_temp + mean_humidity + (1 | year)
##            npar  logLik    AIC    LRT Df Pr(>Chisq)  
## <none>        6 -267.52 547.03                       
## (1 | year)    5 -269.78 549.56 4.5275  1    0.03335 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### So use one model for everything!
### 2022 All together since correlation is low/moderate
a1 <- glm(animals ~ shrub_density + aridity + mean_temp + mean_humidity, data = final_data_2022)
a1
## 
## Call:  glm(formula = animals ~ shrub_density + aridity + mean_temp + 
##     mean_humidity, data = final_data_2022)
## 
## Coefficients:
##   (Intercept)  shrub_density        aridity      mean_temp  mean_humidity  
##       594.701          5.718        -87.742        -15.161         12.542  
## 
## Degrees of Freedom: 23 Total (i.e. Null);  19 Residual
## Null Deviance:       619900 
## Residual Deviance: 339100    AIC: 309.5
anova(a1, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: animals
## 
## Terms added sequentially (first to last)
## 
## 
##               Df Deviance Resid. Df Resid. Dev Pr(>Chi)   
## NULL                             23     619927            
## shrub_density  1    26978        22     592949 0.218915   
## aridity        1    71496        21     521453 0.045349 * 
## mean_temp      1   176488        20     344964 0.001664 **
## mean_humidity  1     5834        19     339131 0.567519   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(a1), "\n")
## Linear Model AIC: 309.4551
#q_model_1 <- glmer(animals ~ shrub_density + aridity  + mean_temp, data = final_data_2022)
#q_model_1
#anova(q_model_1, test = "Chisq")
#cat("Quadratic Model AIC:", AIC(q_model_1), "\n")

if (AIC(quadratic_model_1) < AIC(a1)) {
  cat("Quadratic model is better.\n")
} else {
  cat("Linear model is better or equally good.\n")
}
## Quadratic model is better.
a2 <- glm(richness ~ shrub_density + aridity + mean_temp + mean_humidity, data = final_data_2022)
a2
## 
## Call:  glm(formula = richness ~ shrub_density + aridity + mean_temp + 
##     mean_humidity, data = final_data_2022)
## 
## Coefficients:
##   (Intercept)  shrub_density        aridity      mean_temp  mean_humidity  
##       5.01364        0.05630       -3.01765       -0.08542        0.73972  
## 
## Degrees of Freedom: 23 Total (i.e. Null);  19 Residual
## Null Deviance:       475 
## Residual Deviance: 29.85     AIC: 85.35
anova(a2, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: richness
## 
## Terms added sequentially (first to last)
## 
## 
##               Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                             23     474.96              
## shrub_density  1    6.431        22     468.53 0.0430701 *  
## aridity        1  133.613        21     334.91 < 2.2e-16 ***
## mean_temp      1  284.765        20      50.15 < 2.2e-16 ***
## mean_humidity  1   20.295        19      29.85 0.0003258 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(a2), "\n")
## Linear Model AIC: 85.34806
q_model_2 <- glm(richness ~ shrub_density + I(shrub_density^2) + aridity + I(aridity^2) + mean_temp + I(mean_temp^2) + mean_humidity + I(mean_humidity^2), data = final_data_2022)
q_model_2
## 
## Call:  glm(formula = richness ~ shrub_density + I(shrub_density^2) + 
##     aridity + I(aridity^2) + mean_temp + I(mean_temp^2) + mean_humidity + 
##     I(mean_humidity^2), data = final_data_2022)
## 
## Coefficients:
##        (Intercept)       shrub_density  I(shrub_density^2)             aridity  
##           107.7409             -0.8030              0.0688             99.3440  
##       I(aridity^2)           mean_temp      I(mean_temp^2)       mean_humidity  
##           -10.8127            -37.2859              0.8090              7.0335  
## I(mean_humidity^2)  
##            -0.1173  
## 
## Degrees of Freedom: 23 Total (i.e. Null);  15 Residual
## Null Deviance:       475 
## Residual Deviance: 15.68     AIC: 77.89
anova(q_model_2, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: richness
## 
## Terms added sequentially (first to last)
## 
## 
##                    Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                  23     474.96              
## shrub_density       1    6.431        22     468.53  0.013127 *  
## I(shrub_density^2)  1  106.430        21     362.10 < 2.2e-16 ***
## aridity             1  101.019        20     261.08 < 2.2e-16 ***
## I(aridity^2)        1  232.882        19      28.20 < 2.2e-16 ***
## mean_temp           1    0.134        18      28.06  0.720348    
## I(mean_temp^2)      1   10.796        17      17.27  0.001311 ** 
## mean_humidity       1    0.371        16      16.89  0.551123    
## I(mean_humidity^2)  1    1.215        15      15.68  0.281016    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(q_model_2), "\n")
## Quadratic Model AIC: 77.89283
if (AIC(q_model_2) < AIC(a2)) {
  cat("Quadratic model is better.\n")
} else {
  cat("Linear model is better or equally good.\n")
}
## Quadratic model is better.
a3 <- glm(Average_Evenness ~ shrub_density + aridity + mean_temp + mean_humidity, data = final_data_2022)
a3
## 
## Call:  glm(formula = Average_Evenness ~ shrub_density + aridity + mean_temp + 
##     mean_humidity, data = final_data_2022)
## 
## Coefficients:
##   (Intercept)  shrub_density        aridity      mean_temp  mean_humidity  
##     -0.983958       0.001086      -0.033891       0.026940       0.024481  
## 
## Degrees of Freedom: 23 Total (i.e. Null);  19 Residual
## Null Deviance:       0.1455 
## Residual Deviance: 0.01438   AIC: -97.98
anova(a3, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: Average_Evenness
## 
## Terms added sequentially (first to last)
## 
## 
##               Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                             23   0.145519              
## shrub_density  1 0.007333        22   0.138186  0.001851 ** 
## aridity        1 0.029971        21   0.108215 3.098e-10 ***
## mean_temp      1 0.071611        20   0.036605 < 2.2e-16 ***
## mean_humidity  1 0.022229        19   0.014376 5.950e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(a3), "\n")
## Linear Model AIC: -97.97747
q_model_3 <- glm(Average_Evenness ~ shrub_density + I(shrub_density^2) + aridity + I(aridity^2) + mean_temp + I(mean_temp^2) + mean_humidity + I(mean_humidity^2), data = final_data_2022)
q_model_3
## 
## Call:  glm(formula = Average_Evenness ~ shrub_density + I(shrub_density^2) + 
##     aridity + I(aridity^2) + mean_temp + I(mean_temp^2) + mean_humidity + 
##     I(mean_humidity^2), data = final_data_2022)
## 
## Coefficients:
##        (Intercept)       shrub_density  I(shrub_density^2)             aridity  
##           2.766661           -0.026205            0.002274           -1.027071  
##       I(aridity^2)           mean_temp      I(mean_temp^2)       mean_humidity  
##           0.107674           -0.241692            0.005770            0.146552  
## I(mean_humidity^2)  
##          -0.002244  
## 
## Degrees of Freedom: 23 Total (i.e. Null);  15 Residual
## Null Deviance:       0.1455 
## Residual Deviance: 0.005106  AIC: -114.8
anova(q_model_3, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: Average_Evenness
## 
## Terms added sequentially (first to last)
## 
## 
##                    Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                  23   0.145519              
## shrub_density       1 0.007333        22   0.138186 3.460e-06 ***
## I(shrub_density^2)  1 0.056383        21   0.081804 < 2.2e-16 ***
## aridity             1 0.018828        20   0.062975 1.026e-13 ***
## I(aridity^2)        1 0.049621        19   0.013354 < 2.2e-16 ***
## mean_temp           1 0.001140        18   0.012214  0.067244 .  
## I(mean_temp^2)      1 0.003081        17   0.009133  0.002626 ** 
## mean_humidity       1 0.003583        16   0.005550  0.001176 ** 
## I(mean_humidity^2)  1 0.000444        15   0.005106  0.253210    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(q_model_3), "\n")
## Quadratic Model AIC: -114.822
if (AIC(q_model_3) < AIC(a3)) {
  cat("Quadratic model is better.\n")
} else {
  cat("Linear model is better or equally good.\n")
}
## Quadratic model is better.
b1 <- glm(animals ~ shrub_density + aridity + mean_temp + mean_humidity, data = final_data_2023)
b1
## 
## Call:  glm(formula = animals ~ shrub_density + aridity + mean_temp + 
##     mean_humidity, data = final_data_2023)
## 
## Coefficients:
##   (Intercept)  shrub_density        aridity      mean_temp  mean_humidity  
##      194.4432         1.3948         3.4741        -8.0000        -0.6446  
## 
## Degrees of Freedom: 21 Total (i.e. Null);  17 Residual
## Null Deviance:       12650 
## Residual Deviance: 5738  AIC: 196.8
anova(b1, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: animals
## 
## Terms added sequentially (first to last)
## 
## 
##               Df Deviance Resid. Df Resid. Dev Pr(>Chi)   
## NULL                             21    12648.8            
## shrub_density  1  2257.72        20    10391.1 0.009702 **
## aridity        1  1570.44        19     8820.6 0.031007 * 
## mean_temp      1  3018.31        18     5802.3 0.002787 **
## mean_humidity  1    64.11        17     5738.2 0.662961   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(a1), "\n")
## Linear Model AIC: 309.4551
q_model_1b <- glm(animals ~ shrub_density + I(shrub_density^2) + aridity + I(aridity^2) + mean_temp + I(mean_temp^2) + mean_humidity + I(mean_humidity^2), data = final_data_2023)
q_model_1b
## 
## Call:  glm(formula = animals ~ shrub_density + I(shrub_density^2) + 
##     aridity + I(aridity^2) + mean_temp + I(mean_temp^2) + mean_humidity + 
##     I(mean_humidity^2), data = final_data_2023)
## 
## Coefficients:
##        (Intercept)       shrub_density  I(shrub_density^2)             aridity  
##          -486.3157            -10.0215              0.7863            224.2800  
##       I(aridity^2)           mean_temp      I(mean_temp^2)       mean_humidity  
##           -12.2965            -73.3075              1.4164             30.5245  
## I(mean_humidity^2)  
##            -0.3424  
## 
## Degrees of Freedom: 21 Total (i.e. Null);  13 Residual
## Null Deviance:       12650 
## Residual Deviance: 2388  AIC: 185.5
anova(q_model_1b, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: animals
## 
## Terms added sequentially (first to last)
## 
## 
##                    Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                  21    12648.8              
## shrub_density       1  2257.72        20    10391.1 0.0004549 ***
## I(shrub_density^2)  1    55.46        19    10335.6 0.5826786    
## aridity             1  1622.60        18     8713.0 0.0029563 ** 
## I(aridity^2)        1  1189.72        17     7523.3 0.0109254 *  
## mean_temp           1  2572.77        16     4950.5 0.0001821 ***
## I(mean_temp^2)      1  1451.04        15     3499.5 0.0049430 ** 
## mean_humidity       1   529.64        14     2969.8 0.0894861 .  
## I(mean_humidity^2)  1   582.08        13     2387.7 0.0750425 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(q_model_1b), "\n")
## Quadratic Model AIC: 185.5486
if (AIC(q_model_1b) < AIC(b1)) {
  cat("Quadratic model is better.\n")
} else {
  cat("Linear model is better or equally good.\n")
}
## Quadratic model is better.
b2 <- glm(richness ~ shrub_density + aridity + mean_temp + mean_humidity, data = final_data_2023)
b2
## 
## Call:  glm(formula = richness ~ shrub_density + aridity + mean_temp + 
##     mean_humidity, data = final_data_2023)
## 
## Coefficients:
##   (Intercept)  shrub_density        aridity      mean_temp  mean_humidity  
##      -13.0027         0.1575         0.2489         0.4385         0.1494  
## 
## Degrees of Freedom: 21 Total (i.e. Null);  17 Residual
## Null Deviance:       95.32 
## Residual Deviance: 49.66     AIC: 92.35
anova(b2, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: richness
## 
## Terms added sequentially (first to last)
## 
## 
##               Df Deviance Resid. Df Resid. Dev Pr(>Chi)   
## NULL                             21     95.318            
## shrub_density  1  25.4984        20     69.820 0.003133 **
## aridity        1  16.3249        19     53.495 0.018084 * 
## mean_temp      1   0.3848        18     53.110 0.716651   
## mean_humidity  1   3.4457        17     49.664 0.277465   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(b2), "\n")
## Linear Model AIC: 92.34667
q_model_2b <- glm(richness ~ shrub_density + I(shrub_density^2) + aridity + I(aridity^2) + mean_temp + I(mean_temp^2) + mean_humidity + I(mean_humidity^2), data = final_data_2023)
q_model_2b
## 
## Call:  glm(formula = richness ~ shrub_density + I(shrub_density^2) + 
##     aridity + I(aridity^2) + mean_temp + I(mean_temp^2) + mean_humidity + 
##     I(mean_humidity^2), data = final_data_2023)
## 
## Coefficients:
##        (Intercept)       shrub_density  I(shrub_density^2)             aridity  
##         -172.23925            -0.98555             0.08180            24.28498  
##       I(aridity^2)           mean_temp      I(mean_temp^2)       mean_humidity  
##           -1.33580             2.07678            -0.04053             2.86793  
## I(mean_humidity^2)  
##           -0.02884  
## 
## Degrees of Freedom: 21 Total (i.e. Null);  13 Residual
## Null Deviance:       95.32 
## Residual Deviance: 26.03     AIC: 86.13
anova(q_model_2b, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: richness
## 
## Terms added sequentially (first to last)
## 
## 
##                    Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                  21     95.318              
## shrub_density       1  25.4984        20     69.820 0.0003585 ***
## I(shrub_density^2)  1   0.4031        19     69.417 0.6536249    
## aridity             1  17.3499        18     52.067 0.0032410 ** 
## I(aridity^2)        1   0.0894        17     51.977 0.8326496    
## mean_temp           1   0.3069        16     51.670 0.6953844    
## I(mean_temp^2)      1   8.5804        15     43.090 0.0384262 *  
## mean_humidity       1  12.9355        14     30.155 0.0110236 *  
## I(mean_humidity^2)  1   4.1294        13     26.025 0.1509434    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(q_model_2b), "\n")
## Quadratic Model AIC: 86.12974
if (AIC(q_model_2b) < AIC(b2)) {
  cat("Quadratic model is better.\n")
} else {
  cat("Linear model is better or equally good.\n")
}
## Quadratic model is better.
b3 <- glm(Average_Evenness ~ shrub_density + aridity + mean_temp + mean_humidity, data = final_data_2023)
b3
## 
## Call:  glm(formula = Average_Evenness ~ shrub_density + aridity + mean_temp + 
##     mean_humidity, data = final_data_2023)
## 
## Coefficients:
##   (Intercept)  shrub_density        aridity      mean_temp  mean_humidity  
##      0.612360       0.008000       0.004847      -0.009961      -0.011192  
## 
## Degrees of Freedom: 21 Total (i.e. Null);  17 Residual
## Null Deviance:       0.1608 
## Residual Deviance: 0.01975   AIC: -79.91
anova(b3, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: Average_Evenness
## 
## Terms added sequentially (first to last)
## 
## 
##               Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                             21   0.160797              
## shrub_density  1 0.013715        20   0.147082 0.0005914 ***
## aridity        1 0.039182        19   0.107900 6.369e-09 ***
## mean_temp      1 0.068815        18   0.039085 1.410e-14 ***
## mean_humidity  1 0.019331        17   0.019754 4.530e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(b3), "\n")
## Linear Model AIC: -79.90611
q_model_3b <- glm(Average_Evenness ~ shrub_density + I(shrub_density^2) + aridity + I(aridity^2) + mean_temp + I(mean_temp^2) + mean_humidity + I(mean_humidity^2), data = final_data_2023)
q_model_3b
## 
## Call:  glm(formula = Average_Evenness ~ shrub_density + I(shrub_density^2) + 
##     aridity + I(aridity^2) + mean_temp + I(mean_temp^2) + mean_humidity + 
##     I(mean_humidity^2), data = final_data_2023)
## 
## Coefficients:
##        (Intercept)       shrub_density  I(shrub_density^2)             aridity  
##          6.433e+00           5.868e-03          -1.715e-04          -3.056e-02  
##       I(aridity^2)           mean_temp      I(mean_temp^2)       mean_humidity  
##          1.657e-03          -5.691e-01           1.279e-02          -1.462e-03  
## I(mean_humidity^2)  
##          4.833e-05  
## 
## Degrees of Freedom: 21 Total (i.e. Null);  13 Residual
## Null Deviance:       0.1608 
## Residual Deviance: 0.0002097     AIC: -171.9
anova(q_model_3b, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: Average_Evenness
## 
## Terms added sequentially (first to last)
## 
## 
##                    Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                  21   0.160797              
## shrub_density       1 0.013715        20   0.147082 < 2.2e-16 ***
## I(shrub_density^2)  1 0.020159        19   0.126922 < 2.2e-16 ***
## aridity             1 0.023110        18   0.103812 < 2.2e-16 ***
## I(aridity^2)        1 0.096719        17   0.007093 < 2.2e-16 ***
## mean_temp           1 0.001253        16   0.005840 < 2.2e-16 ***
## I(mean_temp^2)      1 0.005458        15   0.000381 < 2.2e-16 ***
## mean_humidity       1 0.000160        14   0.000221  0.001638 ** 
## I(mean_humidity^2)  1 0.000012        13   0.000210  0.396490    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(q_model_3b), "\n")
## Quadratic Model AIC: -171.9054
if (AIC(q_model_3b) < AIC(b3)) {
  cat("Quadratic model is better.\n")
} else {
  cat("Linear model is better or equally good.\n")
}
## Quadratic model is better.
library(sjPlot)
q_model_1 <- lmer(Average_Evenness ~ shrub_density + aridity  + mean_temp + (1|year), data = final_data) ### Generalize linear mixed effect models
q_model_1
## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: Average_Evenness ~ shrub_density + aridity + mean_temp + (1 |  
##     year)
##    Data: final_data
## REML criterion at convergence: -74.5199
## Random effects:
##  Groups   Name        Std.Dev.
##  year     (Intercept) 0.05073 
##  Residual             0.07346 
## Number of obs: 46, groups:  year, 2
## Fixed Effects:
##   (Intercept)  shrub_density        aridity      mean_temp  
##      0.120540       0.003439       0.011743      -0.004518
anova(q_model_1, test = "Chisq")
## Type III Analysis of Variance Table with Satterthwaite's method
##                 Sum Sq  Mean Sq NumDF  DenDF F value   Pr(>F)   
## shrub_density 0.017489 0.017489     1 41.000  3.2410 0.079177 . 
## aridity       0.048741 0.048741     1 36.939  9.0326 0.004745 **
## mean_temp     0.009448 0.009448     1 41.643  1.7509 0.192988   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("LMER Model AIC:", AIC(q_model_1), "\n")
## LMER Model AIC: -62.51993
tab_model(q_model_1)
  Average_Evenness
Predictors Estimates CI p
(Intercept) 0.12 -0.08 – 0.32 0.228
shrub density 0.00 -0.00 – 0.01 0.079
aridity 0.01 0.00 – 0.02 0.005
mean temp -0.00 -0.01 – 0.00 0.193
Random Effects
σ2 0.01
τ00 year 0.00
ICC 0.32
N year 2
Observations 46
Marginal R2 / Conditional R2 0.276 / 0.510
library(ggplot2)
even_temp <- ggplot(final_data, aes(aridity, Average_Evenness)) +
geom_smooth(method = "lm", formula = y ~ x +I(x^2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") +  labs(x = "aridity", y = "Evenness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
  theme(axis.title.x = element_blank()) + labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
even_temp

cor.test(final_data_2022$aridity, final_data_2022$mean_humidity,method = c("pearson"))
## 
##  Pearson's product-moment correlation
## 
## data:  final_data_2022$aridity and final_data_2022$mean_humidity
## t = 6.3135, df = 22, p-value = 2.36e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5905535 0.9110920
## sample estimates:
##       cor 
## 0.8027203
q_model_2 <- lmer(animals ~ shrub_density + aridity  + mean_temp + (1|year), data = final_data) ### Generalize linear mixed effect models
q_model_2
## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: animals ~ shrub_density + aridity + mean_temp + (1 | year)
##    Data: final_data
## REML criterion at convergence: 535.3932
## Random effects:
##  Groups   Name        Std.Dev.
##  year     (Intercept) 123.4   
##  Residual             103.4   
## Number of obs: 46, groups:  year, 2
## Fixed Effects:
##   (Intercept)  shrub_density        aridity      mean_temp  
##       489.791          3.422          4.482        -18.519
anova(q_model_2, test = "Chisq")
## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## shrub_density  17315   17315     1 41.000   1.621 0.2101201    
## aridity         6761    6761     1 41.958   0.633 0.4307392    
## mean_temp     157818  157818     1 41.275  14.775 0.0004111 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("LMER Model AIC:", AIC(q_model_2), "\n")
## LMER Model AIC: 547.3932
tab_model(q_model_2)
  animals
Predictors Estimates CI p
(Intercept) 489.79 175.18 – 804.41 0.003
shrub density 3.42 -2.01 – 8.85 0.210
aridity 4.48 -6.90 – 15.87 0.431
mean temp -18.52 -28.26 – -8.78 <0.001
Random Effects
σ2 10681.71
τ00 year 15221.82
ICC 0.59
N year 2
Observations 46
Marginal R2 / Conditional R2 0.169 / 0.657
q_model_3 <- lmer(richness ~ shrub_density + aridity  + mean_temp + mean_humidity + (1|year), data = final_data) ### Generalize linear mixed effect models
q_model_3
## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: richness ~ shrub_density + aridity + mean_temp + mean_humidity +  
##     (1 | year)
##    Data: final_data
## REML criterion at convergence: 222.1085
## Random effects:
##  Groups   Name        Std.Dev.
##  year     (Intercept) 2.725   
##  Residual             2.449   
## Number of obs: 46, groups:  year, 2
## Fixed Effects:
##   (Intercept)  shrub_density        aridity      mean_temp  mean_humidity  
##      27.11280        0.12703        0.32017       -0.89528       -0.06545
anova(q_model_3, test = "Chisq")
## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## shrub_density 21.243  21.243     1 40.459  3.5429 0.067006 . 
## aridity       32.049  32.049     1 38.813  5.3452 0.026175 * 
## mean_temp     45.401  45.401     1 40.386  7.5720 0.008835 **
## mean_humidity  1.494   1.494     1 38.798  0.2491 0.620509   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("LMER Model AIC:", AIC(q_model_3), "\n")
## LMER Model AIC: 236.1085
tab_model(q_model_3)
  richness
Predictors Estimates CI p
(Intercept) 27.11 3.21 – 51.01 0.027
shrub density 0.13 -0.01 – 0.26 0.067
aridity 0.32 0.04 – 0.60 0.026
mean temp -0.90 -1.55 – -0.24 0.009
mean humidity -0.07 -0.33 – 0.20 0.620
Random Effects
σ2 6.00
τ00 year 7.43
ICC 0.55
N year 2
Observations 46
Marginal R2 / Conditional R2 0.414 / 0.738
library(emmeans)
arid <- glm(aridity~site * year, data = final_data)
summary(arid)
## 
## Call:
## glm(formula = aridity ~ site * year, data = final_data)
## 
## Coefficients:
##                   Estimate Std. Error    t value Pr(>|t|)    
## (Intercept)     -1.661e+04  2.570e-08 -6.463e+11   <2e-16 ***
## siteCuyama       1.614e+04  3.502e-08  4.608e+11   <2e-16 ***
## siteTecopa      -1.536e+02  3.502e-08 -4.387e+09   <2e-16 ***
## year             8.215e+00  1.271e-11  6.466e+11   <2e-16 ***
## siteCuyama:year -7.980e+00  1.731e-11 -4.609e+11   <2e-16 ***
## siteTecopa:year  7.499e-02  1.731e-11  4.331e+09   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 5.534469e-22)
## 
##     Null deviance: 6.4447e+02  on 45  degrees of freedom
## Residual deviance: 2.2138e-20  on 40  degrees of freedom
## AIC: -2113.4
## 
## Number of Fisher Scoring iterations: 1
anova(arid, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: aridity
## 
## Terms added sequentially (first to last)
## 
## 
##           Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                         45     644.47              
## site       2   137.97        43     506.50 < 2.2e-16 ***
## year       1   339.23        42     167.28 < 2.2e-16 ***
## site:year  2   167.28        40       0.00 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(arid, pairwise~site|year)
## $emmeans
## year = 2022:
##  site    emmean       SE df lower.CL upper.CL
##  Carrizo   5.13 8.32e-12 40     5.13     5.13
##  Cuyama    4.35 8.32e-12 40     4.35     4.35
##  Tecopa    3.12 8.32e-12 40     3.12     3.12
## 
## year = 2023:
##  site    emmean       SE df lower.CL upper.CL
##  Carrizo  13.34 9.60e-12 40    13.34    13.34
##  Cuyama    4.59 8.32e-12 40     4.59     4.59
##  Tecopa   11.41 8.32e-12 40    11.41    11.41
## 
## Confidence level used: 0.95 
## 
## $contrasts
## year = 2022:
##  contrast         estimate       SE df           t.ratio p.value
##  Carrizo - Cuyama    0.778 1.18e-11 40   66120078152.000  <.0001
##  Carrizo - Tecopa    2.010 1.18e-11 40  170900140809.000  <.0001
##  Cuyama - Tecopa     1.232 1.18e-11 40  104780062301.000  <.0001
## 
## year = 2023:
##  contrast         estimate       SE df           t.ratio p.value
##  Carrizo - Cuyama    8.758 1.27e-11 40  689337834957.000  <.0001
##  Carrizo - Tecopa    1.935 1.27e-11 40  152320861822.000  <.0001
##  Cuyama - Tecopa    -6.823 1.18e-11 40 -580044624656.000  <.0001
## 
## P value adjustment: tukey method for comparing a family of 3 estimates
### 
temp <- glm(mean_temp~site * year, data = final_data)
summary(temp)
## 
## Call:
## glm(formula = mean_temp ~ site * year, data = final_data)
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3422.1290   933.9302   3.664 0.000720 ***
## siteCuyama      -2454.3538  1272.7531  -1.928 0.060927 .  
## siteTecopa       6536.5176  1272.7531   5.136 7.68e-06 ***
## year               -1.6806     0.4618  -3.639 0.000774 ***
## siteCuyama:year     1.2132     0.6293   1.928 0.060996 .  
## siteTecopa:year    -3.2291     0.6293  -5.131 7.79e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.7311318)
## 
##     Null deviance: 516.406  on 45  degrees of freedom
## Residual deviance:  29.245  on 40  degrees of freedom
## AIC: 123.71
## 
## Number of Fisher Scoring iterations: 2
anova(temp, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: mean_temp
## 
## Terms added sequentially (first to last)
## 
## 
##           Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                         45     516.41              
## site       2   380.19        43     136.22 < 2.2e-16 ***
## year       1    65.07        42      71.15 < 2.2e-16 ***
## site:year  2    41.91        40      29.25  3.58e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(temp, pairwise~site|year)
## $emmeans
## year = 2022:
##  site    emmean    SE df lower.CL upper.CL
##  Carrizo   24.0 0.302 40     23.4     24.7
##  Cuyama    22.8 0.302 40     22.2     23.4
##  Tecopa    31.4 0.302 40     30.8     32.0
## 
## year = 2023:
##  site    emmean    SE df lower.CL upper.CL
##  Carrizo   22.4 0.349 40     21.7     23.1
##  Cuyama    22.3 0.302 40     21.7     22.9
##  Tecopa    26.5 0.302 40     25.9     27.1
## 
## Confidence level used: 0.95 
## 
## $contrasts
## year = 2022:
##  contrast         estimate    SE df t.ratio p.value
##  Carrizo - Cuyama   1.2539 0.428 40   2.933  0.0150
##  Carrizo - Tecopa  -7.3265 0.428 40 -17.137  <.0001
##  Cuyama - Tecopa   -8.5803 0.428 40 -20.069  <.0001
## 
## year = 2023:
##  contrast         estimate    SE df t.ratio p.value
##  Carrizo - Cuyama   0.0407 0.462 40   0.088  0.9957
##  Carrizo - Tecopa  -4.0974 0.462 40  -8.873  <.0001
##  Cuyama - Tecopa   -4.1380 0.428 40  -9.679  <.0001
## 
## P value adjustment: tukey method for comparing a family of 3 estimates
b1 <- glm(Average_Evenness ~ shrub_density + aridity + mean_temp, data = final_data_2022)
b1
## 
## Call:  glm(formula = Average_Evenness ~ shrub_density + aridity + mean_temp, 
##     data = final_data_2022)
## 
## Coefficients:
##   (Intercept)  shrub_density        aridity      mean_temp  
##      1.074552       0.003864      -0.065522      -0.027385  
## 
## Degrees of Freedom: 23 Total (i.e. Null);  20 Residual
## Null Deviance:       0.1455 
## Residual Deviance: 0.0366    AIC: -77.55
anova(b1, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: Average_Evenness
## 
## Terms added sequentially (first to last)
## 
## 
##               Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                             23   0.145519              
## shrub_density  1 0.007333        22   0.138186   0.04533 *  
## aridity        1 0.029971        21   0.108215 5.195e-05 ***
## mean_temp      1 0.071611        20   0.036605 3.972e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### Use second degree function for paper!
q <- glm(Average_Evenness ~ shrub_density + I(shrub_density^2) + aridity + I(aridity^2) + mean_temp + I(mean_temp^2), data = final_data_2023)
q
## 
## Call:  glm(formula = Average_Evenness ~ shrub_density + I(shrub_density^2) + 
##     aridity + I(aridity^2) + mean_temp + I(mean_temp^2), data = final_data_2023)
## 
## Coefficients:
##        (Intercept)       shrub_density  I(shrub_density^2)             aridity  
##          6.7398056           0.0063338          -0.0001901          -0.0516903  
##       I(aridity^2)           mean_temp      I(mean_temp^2)  
##          0.0028418          -0.5858324           0.0130916  
## 
## Degrees of Freedom: 21 Total (i.e. Null);  15 Residual
## Null Deviance:       0.1608 
## Residual Deviance: 0.0003813     AIC: -162.8
anova(q, test="Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: Average_Evenness
## 
## Terms added sequentially (first to last)
## 
## 
##                    Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                  21   0.160797              
## shrub_density       1 0.013715        20   0.147082 < 2.2e-16 ***
## I(shrub_density^2)  1 0.020159        19   0.126922 < 2.2e-16 ***
## aridity             1 0.023110        18   0.103812 < 2.2e-16 ***
## I(aridity^2)        1 0.096719        17   0.007093 < 2.2e-16 ***
## mean_temp           1 0.001253        16   0.005840 2.191e-12 ***
## I(mean_temp^2)      1 0.005458        15   0.000381 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
corelation_fig <- ggplot(final_data, aes(x = mean_temp, y = mean_humidity)) +
  geom_point() + labs(fill = "", x = "Temperature (°C)", y = "Humidity (%)") +theme_classic()+theme(aspect.ratio = 1) + geom_smooth(method = lm, se = F) + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7)
cor <- corelation_fig + stat_cor(method = "pearson")
cor

### Inflection stuff
library(inflection)
q <- glm(animals ~ shrub_density + I(shrub_density^2) + aridity + I(aridity^2) + mean_temp + I(mean_temp^2), data = final_data_2023)
q
## 
## Call:  glm(formula = animals ~ shrub_density + I(shrub_density^2) + 
##     aridity + I(aridity^2) + mean_temp + I(mean_temp^2), data = final_data_2023)
## 
## Coefficients:
##        (Intercept)       shrub_density  I(shrub_density^2)             aridity  
##          3942.1978             -2.6344              0.2613           -114.4117  
##       I(aridity^2)           mean_temp      I(mean_temp^2)  
##             6.4616           -309.4103              6.7500  
## 
## Degrees of Freedom: 21 Total (i.e. Null);  15 Residual
## Null Deviance:       12650 
## Residual Deviance: 3499  AIC: 190
anova(q, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: animals
## 
## Terms added sequentially (first to last)
## 
## 
##                    Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                  21    12648.8              
## shrub_density       1  2257.72        20    10391.1 0.0018654 ** 
## I(shrub_density^2)  1    55.46        19    10335.6 0.6258710    
## aridity             1  1622.60        18     8713.0 0.0083581 ** 
## I(aridity^2)        1  1189.72        17     7523.3 0.0239314 *  
## mean_temp           1  2572.77        16     4950.5 0.0008975 ***
## I(mean_temp^2)      1  1451.04        15     3499.5 0.0126336 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#abline(q, col = "red")
                  
#inflection_point <- find_inflection(final_data$animals, predict(q))
#cat("Inflection Point:", inflection_point, "\n")



#inflection_point <- find_inflection(your_data$abundance, predict(model))
#cat("Inflection Point:", inflection_point, "\n")
 
abund_dens_2022 <- ggplot(final_data_2022, aes(shrub_density, animals)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") +  labs(x = "Shrub Density 20m Radius", y = "Animal Abundance") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))

### 2022 Density
library(inflection)
f=function(x){(1/8*x+1/2*x^2)/(1+1/8*x+1/2*x^2)}
x = final_data_2022$shrub_density
y = final_data_2022$animals

cc=check_curve(x,y)
cc
## $ctype
## [1] "concave"
## 
## $index
## [1] 1
ese(x,y,cc$index)
##     j1 j2 chi
## ESE 13  8 NaN
ipbese=bese(x,y,cc$index)
ipbese$iplast
## [1] NaN
inflection_point_dens_2022 <- final_data_2022[which.min(final_data_2022$y), ]
Fig1_a <- abund_dens_2022 + geom_vline(xintercept = inflection_point_dens_2022$shrub_density, linetype = 2)
Fig1_a

rich_dens_2022 <- ggplot(final_data_2022, aes(shrub_density, richness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") +  labs(x = "Shrub Density 20m Radius", y = "Richness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
  theme(axis.title.x = element_blank()) + labs(tag = "B") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
# Assuming your data frame is named final_data_2022
# Assuming your data has columns 'shrub_density', 'animals', 'richness', 'evenness'

# Load necessary libraries
# Install and load the bcp package
#install.packages("changepoint")
library(changepoint)

# Function to find inflection point using bcp package
find_inflection_point_cpt <- function(x, y) {
  data_ts <- data.frame(y = y)
  cpt_result <- cpt.var(data_ts$y, method = "BinSeg")
  change_point <- cpts(cpt_result)[1]
  return(change_point)
}

# Fit models and find inflection points using bcp
inflection_animals <- find_inflection_point_cpt(final_data_2022$shrub_density, final_data_2022$animals)
inflection_richness <- find_inflection_point_cpt(final_data_2022$shrub_density, final_data_2022$richness)
inflection_evenness <- find_inflection_point_cpt(final_data_2022$shrub_density, final_data_2022$Average_Evenness)


# Print inflection points
cat("Inflection Point for Animals:", inflection_animals, "\n")
## Inflection Point for Animals: 12
cat("Inflection Point for Richness:", inflection_richness, "\n")
## Inflection Point for Richness: 8
cat("Inflection Point for Evenness:", inflection_evenness, "\n")
## Inflection Point for Evenness: NA
# Plot individual graphs
# ... (ggplot code and plotting)
abund_dens_2022 <- ggplot(final_data_2022, aes(shrub_density, animals)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") +  labs(x = "Shrub Density per20m Radius", y = "Abundance") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
Fig1_a <- abund_dens_2022 + geom_vline(xintercept = inflection_animals, linetype = 2)
Fig1_a

rich_dens_2022 <- ggplot(final_data_2022, aes(shrub_density, richness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") +  labs(x = "Shrub Density per 20m Radius", y = "Richness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
labs(tag = "B") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))

Fig1_b <- rich_dens_2022 + geom_vline(xintercept = inflection_richness, linetype = 2)
Fig1_b

even_dens_2022 <- ggplot(final_data_2022, aes(shrub_density, Average_Evenness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") +  labs(x = "Shrub Density per 20m Radius", y = "Evenness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+
   labs(tag = "C") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
Fig1_c <- even_dens_2022 + geom_vline(xintercept = inflection_evenness, linetype = 2)
Fig1_c

### Now for 2023

find_inflection_point_cpt <- function(x, y) {
  data_ts <- data.frame(y = y)
  cpt_result <- cpt.var(data_ts$y, method = "BinSeg")
  change_point <- cpts(cpt_result)
  return(change_point)
}

# Standardize variables
final_data_2022$standardized_aridity <- scale(final_data_2022$aridity)

# Adjust your model accordingly
inflection_animals <- find_inflection_point_cpt(final_data_2022$standardized_aridity, final_data_2022$animals)
inflection_richness <- find_inflection_point_cpt(final_data_2022$standardized_aridity, final_data_2022$richness)
inflection_evenness <- find_inflection_point_cpt(final_data_2022$standardized_aridity, final_data_2022$Average_Evenness)


# Fit models and find inflection points using changepoint
inflection_animals <- find_inflection_point_cpt(final_data_2022$aridity, final_data_2022$animals)
inflection_richness <- find_inflection_point_cpt(final_data_2022$aridity, final_data_2022$richness)
inflection_evenness <- find_inflection_point_cpt(final_data_2022$aridity, final_data_2022$Average_Evenness)

# Print inflection points
cat("Inflection Point for Animals:", ifelse(length(inflection_animals) > 0, inflection_animals, "No clear inflection point"), "\n")
## Inflection Point for Animals: 12
cat("Inflection Point for Richness:", ifelse(length(inflection_richness) > 0, inflection_richness, "No clear inflection point"), "\n")
## Inflection Point for Richness: 8
cat("Inflection Point for Evenness:", ifelse(length(inflection_evenness) > 0, inflection_evenness, "No clear inflection point"), "\n")
## Inflection Point for Evenness: No clear inflection point
# Function to find inflection point using quadratic model
find_inflection_point_quadratic <- function(x, y) {
  quadratic_model <- lm(y ~ poly(x, 2, raw = TRUE))
  vertex <- -coef(quadratic_model)[2] / (2 * coef(quadratic_model)[3])
  return(vertex)
}

# Find inflection points using quadratic model
den_inflection_animals <- find_inflection_point_quadratic(final_data_2022$shrub_density, final_data_2022$animals)
den_inflection_richness <- find_inflection_point_quadratic(final_data_2022$shrub_density, final_data_2022$richness)
den_inflection_evenness <- find_inflection_point_quadratic(final_data_2022$shrub_density, final_data_2022$Average_Evenness)

# Print inflection points
cat("Inflection Point for Animals:", den_inflection_animals, "\n")
## Inflection Point for Animals: 8.518166
cat("Inflection Point for Richness:", den_inflection_richness, "\n")
## Inflection Point for Richness: 5.595659
cat("Inflection Point for Evenness:", den_inflection_evenness, "\n")
## Inflection Point for Evenness: 5.493175
abund_dens_2022 <- ggplot(final_data_2022, aes(shrub_density, animals)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") +  labs(x = "Shrub Density per 20m Radius", y = "Animal Abundance") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
Fig1_a <- abund_dens_2022 + geom_vline(xintercept = den_inflection_animals, linetype = 2)
Fig1_a

rich_dens_2022 <- ggplot(final_data_2022, aes(shrub_density, richness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") +  labs(x = "Shrub Density per 20m Radius", y = "Richness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
labs(tag = "B") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))

Fig1_b <- rich_dens_2022 + geom_vline(xintercept = den_inflection_richness, linetype = 2)
Fig1_b

even_dens_2022 <- ggplot(final_data_2022, aes(shrub_density, Average_Evenness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") +  labs(x = "Shrub Density per 20m Radius", y = "Evenness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+
   labs(tag = "C") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
Fig1_c <- even_dens_2022 + geom_vline(xintercept = den_inflection_evenness, linetype = 2)
Fig1_c

### Aridity
find_inflection_point_quadratic <- function(x, y) {
  quadratic_model <- lm(y ~ poly(x, 2, raw = TRUE))
  vertex <- -coef(quadratic_model)[2] / (2 * coef(quadratic_model)[3])
  return(vertex)
}

# Find inflection points using quadratic model
arid_inflection_animals <- find_inflection_point_quadratic(final_data_2022$aridity, final_data_2022$animals)
arid_inflection_richness <- find_inflection_point_quadratic(final_data_2022$aridity, final_data_2022$richness)
arid_inflection_evenness <- find_inflection_point_quadratic(final_data_2022$aridity, final_data_2022$Average_Evenness)

# Print inflection points
cat("Inflection Point for Animals:", inflection_animals, "\n")
## Inflection Point for Animals: 12
cat("Inflection Point for Richness:", inflection_richness, "\n")
## Inflection Point for Richness: 8
cat("Inflection Point for Evenness:", inflection_evenness, "\n")
## Inflection Point for Evenness:
abund_arid_2022 <- ggplot(final_data_2022, aes(aridity, animals)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") +  labs(x = "Aridity", y = "Animal Abundance") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))

Fig2_a <- abund_arid_2022 + geom_vline(xintercept = arid_inflection_animals, linetype = 2)
Fig2_a

rich_arid_2022 <- ggplot(final_data_2022, aes(aridity, richness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") +  labs(x = "Aridity", y = "Richness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "B") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))

Fig2_b <- rich_arid_2022 + geom_vline(xintercept = arid_inflection_richness, linetype = 2)
Fig2_b

even_arid_2022 <- ggplot(final_data_2022, aes(aridity, Average_Evenness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") +  labs(x = "Aridity", y = "Evenness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "C") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))

Fig2_c <- even_arid_2022 + geom_vline(xintercept = arid_inflection_evenness, linetype = 2)
Fig2_c

# Temp
find_inflection_point_quadratic <- function(x, y) {
  quadratic_model <- lm(y ~ poly(x, 2, raw = TRUE))
  vertex <- -coef(quadratic_model)[2] / (2 * coef(quadratic_model)[3])
  return(vertex)
}

# Find inflection points using quadratic model
temp_inflection_animals_2022 <- find_inflection_point_quadratic(final_data_2022$mean_temp, final_data_2022$animals)
temp_inflection_richness_2022 <- find_inflection_point_quadratic(final_data_2022$mean_temp, final_data_2022$richness)
temp_inflection_evenness_2022 <- find_inflection_point_quadratic(final_data_2022$mean_temp, final_data_2022$Average_Evenness)

# Print inflection points
cat("Inflection Point for Animals:", temp_inflection_animals_2022, "\n")
## Inflection Point for Animals: 3.261247
cat("Inflection Point for Richness:", temp_inflection_richness_2022, "\n")
## Inflection Point for Richness: 31.81689
cat("Inflection Point for Evenness:", temp_inflection_evenness_2022, "\n")
## Inflection Point for Evenness: 46.82795
### Temperature might be linear not second degree?

abund_temp_2022 <- ggplot(final_data_2022, aes(mean_temp, animals)) +
geom_smooth(method = "lm", se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") +  labs(x = "Temperature (°C)", y = "Abundance") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
abund_temp_2022

Fig3_a <- abund_temp_2022

rich_temp_2022 <- ggplot(final_data_2022, aes(mean_temp, richness)) +
geom_smooth(method = "lm", se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") +  labs(x = "Temperature (°C)", y = "Richness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "B") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
rich_temp_2022

Fig3_b <- rich_temp_2022 

even_temp_2022 <- ggplot(final_data_2022, aes(mean_temp, Average_Evenness)) +
geom_smooth(method = "lm", se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") +  labs(x = "Temperature (°C)", y = "Evenness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "C") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
even_temp_2022

Fig3_c <- even_temp_2022
find_inflection_point_quadratic <- function(x, y) {
  quadratic_model <- lm(y ~ poly(x, 2, raw = TRUE))
  vertex <- -coef(quadratic_model)[2] / (2 * coef(quadratic_model)[3])
  return(vertex)
}

# Find inflection points using quadratic model
den_inflection_animals_2023 <- find_inflection_point_quadratic(final_data_2023$shrub_density, final_data_2023$animals)
den_inflection_richness_2023 <- find_inflection_point_quadratic(final_data_2023$shrub_density, final_data_2023$richness)
den_inflection_evenness_2023 <- find_inflection_point_quadratic(final_data_2023$shrub_density, final_data_2023$Average_Evenness)

# Print inflection points
cat("Inflection Point for Animals:", den_inflection_animals_2023, "\n")
## Inflection Point for Animals: 11.78866
cat("Inflection Point for Richness:", den_inflection_richness_2023, "\n")
## Inflection Point for Richness: 13.25782
cat("Inflection Point for Evenness:", den_inflection_evenness_2023, "\n")
## Inflection Point for Evenness: 6.598329
abund_dens_2023 <- ggplot(final_data_2023, aes(shrub_density, animals)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") +  labs(x = "Shrub Density per 20m Radius", y = "Abundance") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "D") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
Fig1_d <- abund_dens_2023 + geom_vline(xintercept = den_inflection_animals_2023, linetype = 2)
Fig1_d

rich_dens_2023 <- ggplot(final_data_2023, aes(shrub_density, richness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") +  labs(x = "Shrub Density per 20m Radius", y = "Richness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
labs(tag = "E") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))

Fig1_e <- rich_dens_2023 + geom_vline(xintercept = den_inflection_richness_2023, linetype = 2)
Fig1_e

even_dens_2023 <- ggplot(final_data_2023, aes(shrub_density, Average_Evenness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") +  labs(x = "Shrub Density per 20m Radius", y = "Evenness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+
   labs(tag = "F") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
Fig1_f <- even_dens_2023 + geom_vline(xintercept = den_inflection_evenness_2023, linetype = 2)
Fig1_f

### Aridity
find_inflection_point_quadratic <- function(x, y) {
  quadratic_model <- lm(y ~ poly(x, 2, raw = TRUE))
  vertex <- -coef(quadratic_model)[2] / (2 * coef(quadratic_model)[3])
  return(vertex)
}

# Find inflection points using quadratic model
arid_inflection_animals_2023 <- find_inflection_point_quadratic(final_data_2023$aridity, final_data_2023$animals)
arid_inflection_richness_2023 <- find_inflection_point_quadratic(final_data_2023$aridity, final_data_2023$richness)
arid_inflection_evenness_2023 <- find_inflection_point_quadratic(final_data_2023$aridity, final_data_2023$Average_Evenness)

# Print inflection points
cat("Inflection Point for Animals:", arid_inflection_animals_2023, "\n")
## Inflection Point for Animals: 7.870704
cat("Inflection Point for Richness:", arid_inflection_richness_2023, "\n")
## Inflection Point for Richness: 5.866049
cat("Inflection Point for Evenness:", arid_inflection_evenness_2023, "\n")
## Inflection Point for Evenness: 9.157136
abund_arid_2022 <- ggplot(final_data_2023, aes(aridity, animals)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") +  labs(x = "Aridity", y = "Animal Abundance") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))

Fig2_d <- abund_arid_2022 + geom_vline(xintercept = arid_inflection_animals_2023, linetype = 2)
Fig2_d

rich_arid_2022 <- ggplot(final_data_2023, aes(aridity, richness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") +  labs(x = "Aridity", y = "Richness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))

Fig2_e <- rich_arid_2022 + geom_vline(xintercept = arid_inflection_richness_2023, linetype = 2)
Fig2_e

even_arid_2022 <- ggplot(final_data_2023, aes(aridity, Average_Evenness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") +  labs(x = "Aridity", y = "Evenness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))

Fig2_f <- even_arid_2022 + geom_vline(xintercept = arid_inflection_evenness_2023, linetype = 2)
Fig2_f

### Temp
find_inflection_point_quadratic <- function(x, y) {
  quadratic_model <- lm(y ~ poly(x, 2, raw = TRUE))
  vertex <- -coef(quadratic_model)[2] / (2 * coef(quadratic_model)[3])
  return(vertex)
}

# Find inflection points using quadratic model
temp_inflection_animals_2023 <- find_inflection_point_quadratic(final_data_2023$mean_temp, final_data_2023$animals)
temp_inflection_richness_2023 <- find_inflection_point_quadratic(final_data_2023$mean_temp, final_data_2023$richness)
temp_inflection_evenness_2023 <- find_inflection_point_quadratic(final_data_2023$mean_temp, final_data_2023$Average_Evenness)

# Print inflection points
cat("Inflection Point for Animals:", temp_inflection_animals_2023, "\n")
## Inflection Point for Animals: 24.45404
cat("Inflection Point for Richness:", temp_inflection_richness_2023, "\n")
## Inflection Point for Richness: 23.81655
cat("Inflection Point for Evenness:", temp_inflection_evenness_2023, "\n")
## Inflection Point for Evenness: 22.71232
abund_temp_2023 <- ggplot(final_data_2023, aes(mean_temp, animals)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") +  labs(x = "Temperature (°C)", y = "Abundance") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "D") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))

Fig3_d <- abund_temp_2023 + geom_vline(xintercept = temp_inflection_animals_2023, linetype = 2)
Fig3_d

rich_temp_2023 <- ggplot(final_data_2023, aes(mean_temp, richness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") +  labs(x = "Temperature (°C)", y = "Richness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "E") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))

Fig3_e <- rich_temp_2023 + geom_vline(xintercept = temp_inflection_richness_2023, linetype = 2)
Fig3_e

even_temp_2023 <- ggplot(final_data_2023, aes(mean_temp, Average_Evenness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") +  labs(x = "Temperature (°C)", y = "Evenness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "F") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
even_temp_2023

Fig3_f <- even_temp_2023 + geom_vline(xintercept = temp_inflection_evenness_2023, linetype = 2)
Fig3_f

library(patchwork)
Figure_1 <- (Fig1_a + Fig1_b + Fig1_c) / (Fig1_d + Fig1_e + Fig1_f)
Figure_1

combined_plots <- (
  Fig1_a + Fig1_d + Fig1_b + Fig1_e + Fig1_c + Fig1_f
) + plot_layout(ncol = 2)
combined_plots

combined_plots_temp <- (Fig3_a + Fig3_d + Fig3_b + Fig3_e + Fig3_c +Fig3_f) + plot_layout(ncol=2)
combined_plots_temp

m1.1 <- glm(animals ~ poly(shrub_density,2, raw = TRUE) + mean_temp, family = poisson, final_data_2022)
summary(m1.1)
## 
## Call:
## glm(formula = animals ~ poly(shrub_density, 2, raw = TRUE) + 
##     mean_temp, family = poisson, data = final_data_2022)
## 
## Coefficients:
##                                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                         15.60486    0.34063   45.81   <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)1  0.57397    0.02423   23.68   <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)2 -0.04175    0.00203  -20.57   <2e-16 ***
## mean_temp                           -0.46818    0.01508  -31.05   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 4036.2  on 23  degrees of freedom
## Residual deviance: 1108.0  on 20  degrees of freedom
## AIC: 1243.9
## 
## Number of Fisher Scoring iterations: 5
a <- anova(m1.1, test = "Chisq")
a
## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: animals
## 
## Terms added sequentially (first to last)
## 
## 
##                                    Df Deviance Resid. Df Resid. Dev  Pr(>Chi)
## NULL                                                  23     4036.2          
## poly(shrub_density, 2, raw = TRUE)  2   227.94        21     3808.2 < 2.2e-16
## mean_temp                           1  2700.17        20     1108.0 < 2.2e-16
##                                       
## NULL                                  
## poly(shrub_density, 2, raw = TRUE) ***
## mean_temp                          ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1.2 <- glm(animals ~  poly(shrub_density, 2, raw = TRUE) + mean_temp, family = poisson, final_data_2023)
summary(m1.2)
## 
## Call:
## glm(formula = animals ~ poly(shrub_density, 2, raw = TRUE) + 
##     mean_temp, family = poisson, data = final_data_2023)
## 
## Coefficients:
##                                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                          6.111428   0.482175  12.675  < 2e-16 ***
## poly(shrub_density, 2, raw = TRUE)1  0.217399   0.046585   4.667 3.06e-06 ***
## poly(shrub_density, 2, raw = TRUE)2 -0.012718   0.003905  -3.257  0.00113 ** 
## mean_temp                           -0.147569   0.019964  -7.392 1.45e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 468.53  on 21  degrees of freedom
## Residual deviance: 311.28  on 18  degrees of freedom
## AIC: 416.32
## 
## Number of Fisher Scoring iterations: 5
anova(m1.2, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: animals
## 
## Terms added sequentially (first to last)
## 
## 
##                                    Df Deviance Resid. Df Resid. Dev  Pr(>Chi)
## NULL                                                  21     468.53          
## poly(shrub_density, 2, raw = TRUE)  2  101.012        19     367.52 < 2.2e-16
## mean_temp                           1   56.245        18     311.28 6.399e-14
##                                       
## NULL                                  
## poly(shrub_density, 2, raw = TRUE) ***
## mean_temp                          ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m2.1 <- glm(richness ~  poly(shrub_density, 2, raw = TRUE) + mean_temp, family = gaussian, final_data_2022)
summary(m2.1)
## 
## Call:
## glm(formula = richness ~ poly(shrub_density, 2, raw = TRUE) + 
##     mean_temp, family = gaussian, data = final_data_2022)
## 
## Coefficients:
##                                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                         29.58881    3.29609   8.977 1.88e-08 ***
## poly(shrub_density, 2, raw = TRUE)1 -1.02524    0.55688  -1.841   0.0805 .  
## poly(shrub_density, 2, raw = TRUE)2  0.09754    0.04732   2.061   0.0525 .  
## mean_temp                           -0.91159    0.12537  -7.271 4.93e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 4.968847)
## 
##     Null deviance: 474.958  on 23  degrees of freedom
## Residual deviance:  99.377  on 20  degrees of freedom
## AIC: 112.21
## 
## Number of Fisher Scoring iterations: 2
anova(m2.1, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: richness
## 
## Terms added sequentially (first to last)
## 
## 
##                                    Df Deviance Resid. Df Resid. Dev  Pr(>Chi)
## NULL                                                  23     474.96          
## poly(shrub_density, 2, raw = TRUE)  2   112.86        21     362.10 1.169e-05
## mean_temp                           1   262.72        20      99.38 3.557e-13
##                                       
## NULL                                  
## poly(shrub_density, 2, raw = TRUE) ***
## mean_temp                          ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m2.2 <- glm(richness ~ poly(shrub_density, 2, raw = TRUE) + mean_temp, family = gaussian, final_data_2023)
summary(m2.2)
## 
## Call:
## glm(formula = richness ~ poly(shrub_density, 2, raw = TRUE) + 
##     mean_temp, family = gaussian, data = final_data_2023)
## 
## Coefficients:
##                                      Estimate Std. Error t value Pr(>|t|)
## (Intercept)                          2.434916   4.820063   0.505    0.620
## poly(shrub_density, 2, raw = TRUE)1  0.300793   0.473720   0.635    0.533
## poly(shrub_density, 2, raw = TRUE)2 -0.009007   0.040410  -0.223    0.826
## mean_temp                            0.081306   0.197500   0.412    0.685
## 
## (Dispersion parameter for gaussian family taken to be 3.820508)
## 
##     Null deviance: 95.318  on 21  degrees of freedom
## Residual deviance: 68.769  on 18  degrees of freedom
## AIC: 97.507
## 
## Number of Fisher Scoring iterations: 2
anova(m2.2, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: richness
## 
## Terms added sequentially (first to last)
## 
## 
##                                    Df Deviance Resid. Df Resid. Dev Pr(>Chi)  
## NULL                                                  21     95.318           
## poly(shrub_density, 2, raw = TRUE)  2  25.9015        19     69.417  0.03372 *
## mean_temp                           1   0.6475        18     68.769  0.68058  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m3.1 <- glm(Average_Evenness ~  poly(shrub_density, 2, raw = TRUE) + mean_temp, family = gaussian, final_data_2022)
summary(m3.1)
## 
## Call:
## glm(formula = Average_Evenness ~ poly(shrub_density, 2, raw = TRUE) + 
##     mean_temp, family = gaussian, data = final_data_2022)
## 
## Coefficients:
##                                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                          0.4201571  0.0592359   7.093  7.1e-07 ***
## poly(shrub_density, 2, raw = TRUE)1 -0.0346305  0.0100080  -3.460  0.00247 ** 
## poly(shrub_density, 2, raw = TRUE)2  0.0032634  0.0008504   3.838  0.00103 ** 
## mean_temp                           -0.0125390  0.0022530  -5.565  1.9e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.001604821)
## 
##     Null deviance: 0.145519  on 23  degrees of freedom
## Residual deviance: 0.032096  on 20  degrees of freedom
## AIC: -80.701
## 
## Number of Fisher Scoring iterations: 2
anova(m3.1, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: Average_Evenness
## 
## Terms added sequentially (first to last)
## 
## 
##                                    Df Deviance Resid. Df Resid. Dev  Pr(>Chi)
## NULL                                                  23   0.145519          
## poly(shrub_density, 2, raw = TRUE)  2 0.063715        21   0.081804 2.392e-09
## mean_temp                           1 0.049707        20   0.032096 2.616e-08
##                                       
## NULL                                  
## poly(shrub_density, 2, raw = TRUE) ***
## mean_temp                          ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m3.2 <- glm(Average_Evenness ~ poly(shrub_density, 2, raw = TRUE) + mean_temp, family = gaussian, final_data_2023)
summary(m3.2)
## 
## Call:
## glm(formula = Average_Evenness ~ poly(shrub_density, 2, raw = TRUE) + 
##     mean_temp, family = gaussian, data = final_data_2023)
## 
## Coefficients:
##                                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                         -0.6385620  0.1181717  -5.404 3.91e-05 ***
## poly(shrub_density, 2, raw = TRUE)1  0.0238482  0.0116140   2.053   0.0549 .  
## poly(shrub_density, 2, raw = TRUE)2 -0.0014857  0.0009907  -1.500   0.1510    
## mean_temp                            0.0295605  0.0048420   6.105 9.10e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.002296376)
## 
##     Null deviance: 0.160797  on 21  degrees of freedom
## Residual deviance: 0.041335  on 18  degrees of freedom
## AIC: -65.663
## 
## Number of Fisher Scoring iterations: 2
anova(m3.2, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: Average_Evenness
## 
## Terms added sequentially (first to last)
## 
## 
##                                    Df Deviance Resid. Df Resid. Dev  Pr(>Chi)
## NULL                                                  21   0.160797          
## poly(shrub_density, 2, raw = TRUE)  2 0.033874        19   0.126922 0.0006263
## mean_temp                           1 0.085588        18   0.041335 1.028e-09
##                                       
## NULL                                  
## poly(shrub_density, 2, raw = TRUE) ***
## mean_temp                          ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(tidyverse)
final_data <- read.csv("Final Data.csv")


site <- final_data %>%
  dplyr::select(site_code, year, shrub_density, animals, richness, Average_Evenness, mean_temp, aridity) %>%
  group_by(site_code, year) %>%
  summarise(shrub_density = sum(shrub_density), animals = mean(animals), richness = mean(richness), evenness = mean(Average_Evenness), aridity = aridity)%>%
    mutate(animals = as.integer(animals))

site_2022 <- site%>%
  filter(year == 2022)

site_2023 <- site %>%
  filter(year ==2023)
m4 <- glm(animals ~  poly(shrub_density,2, raw = TRUE) * aridity, family = poisson, site_2022)
summary(m4)
## 
## Call:
## glm(formula = animals ~ poly(shrub_density, 2, raw = TRUE) * 
##     aridity, family = poisson, data = site_2022)
## 
## Coefficients:
##                                               Estimate Std. Error z value
## (Intercept)                                  7.423e+00  5.807e+04   0.000
## poly(shrub_density, 2, raw = TRUE)1         -9.742e-01  2.046e+03   0.000
## poly(shrub_density, 2, raw = TRUE)2          5.529e-02  2.561e+01   0.002
## aridity                                     -1.018e+01  2.074e+02  -0.049
## poly(shrub_density, 2, raw = TRUE)1:aridity  1.074e+00  6.136e+02   0.002
## poly(shrub_density, 2, raw = TRUE)2:aridity -3.099e-02  2.801e+01  -0.001
##                                             Pr(>|z|)
## (Intercept)                                    1.000
## poly(shrub_density, 2, raw = TRUE)1            1.000
## poly(shrub_density, 2, raw = TRUE)2            0.998
## aridity                                        0.961
## poly(shrub_density, 2, raw = TRUE)1:aridity    0.999
## poly(shrub_density, 2, raw = TRUE)2:aridity    0.999
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 3.2932e+03  on 23  degrees of freedom
## Residual deviance: 2.2327e-10  on 18  degrees of freedom
## AIC: 140.45
## 
## Number of Fisher Scoring iterations: 22
anova(m4, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: animals
## 
## Terms added sequentially (first to last)
## 
## 
##                                            Df Deviance Resid. Df Resid. Dev
## NULL                                                          23     3293.2
## poly(shrub_density, 2, raw = TRUE)          2  2223.43        21     1069.8
## aridity                                     1   884.14        20      185.6
## poly(shrub_density, 2, raw = TRUE):aridity  2   185.64        18        0.0
##                                             Pr(>Chi)    
## NULL                                                    
## poly(shrub_density, 2, raw = TRUE)         < 2.2e-16 ***
## aridity                                    < 2.2e-16 ***
## poly(shrub_density, 2, raw = TRUE):aridity < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m5 <- glm(richness ~  poly(shrub_density,2, raw = TRUE) * aridity, family = gaussian, site_2022)
summary(m5)
## 
## Call:
## glm(formula = richness ~ poly(shrub_density, 2, raw = TRUE) * 
##     aridity, family = gaussian, data = site_2022)
## 
## Coefficients:
##                                               Estimate Std. Error    t value
## (Intercept)                                 -2.903e+02  1.642e-11 -1.768e+13
## poly(shrub_density, 2, raw = TRUE)1          2.881e+01  1.373e-12  2.098e+13
## poly(shrub_density, 2, raw = TRUE)2         -6.272e-01  2.876e-14 -2.181e+13
## aridity                                      9.336e+01  5.266e-12  1.773e+13
## poly(shrub_density, 2, raw = TRUE)1:aridity -8.689e+00  4.420e-13 -1.966e+13
## poly(shrub_density, 2, raw = TRUE)2:aridity  1.882e-01  9.266e-15  2.032e+13
##                                             Pr(>|t|)    
## (Intercept)                                   <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)1           <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)2           <2e-16 ***
## aridity                                       <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)1:aridity   <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)2:aridity   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 5.024603e-26)
## 
##     Null deviance: 4.5921e+02  on 23  degrees of freedom
## Residual deviance: 9.0443e-25  on 18  degrees of freedom
## AIC: -1322.9
## 
## Number of Fisher Scoring iterations: 1
anova(m5, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: richness
## 
## Terms added sequentially (first to last)
## 
## 
##                                            Df Deviance Resid. Df Resid. Dev
## NULL                                                          23     459.21
## poly(shrub_density, 2, raw = TRUE)          2  293.998        21     165.21
## aridity                                     1  134.733        20      30.48
## poly(shrub_density, 2, raw = TRUE):aridity  2   30.477        18       0.00
##                                             Pr(>Chi)    
## NULL                                                    
## poly(shrub_density, 2, raw = TRUE)         < 2.2e-16 ***
## aridity                                    < 2.2e-16 ***
## poly(shrub_density, 2, raw = TRUE):aridity < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m6 <- glm(evenness ~  poly(shrub_density,2, raw = TRUE) * aridity, family = gaussian, site_2022)
summary(m6)
## 
## Call:
## glm(formula = evenness ~ poly(shrub_density, 2, raw = TRUE) * 
##     aridity, family = gaussian, data = site_2022)
## 
## Coefficients:
##                                               Estimate Std. Error    t value
## (Intercept)                                 -8.772e+00  4.329e-13 -2.026e+13
## poly(shrub_density, 2, raw = TRUE)1          8.151e-01  3.621e-14  2.251e+13
## poly(shrub_density, 2, raw = TRUE)2         -1.775e-02  7.583e-16 -2.341e+13
## aridity                                      2.823e+00  1.388e-13  2.034e+13
## poly(shrub_density, 2, raw = TRUE)1:aridity -2.539e-01  1.165e-14 -2.179e+13
## poly(shrub_density, 2, raw = TRUE)2:aridity  5.514e-03  2.443e-16  2.257e+13
##                                             Pr(>|t|)    
## (Intercept)                                   <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)1           <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)2           <2e-16 ***
## aridity                                       <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)1:aridity   <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)2:aridity   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 3.493112e-29)
## 
##     Null deviance: 1.3525e-01  on 23  degrees of freedom
## Residual deviance: 6.2876e-28  on 18  degrees of freedom
## AIC: -1497.4
## 
## Number of Fisher Scoring iterations: 1
anova(m6, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: evenness
## 
## Terms added sequentially (first to last)
## 
## 
##                                            Df Deviance Resid. Df Resid. Dev
## NULL                                                          23   0.135252
## poly(shrub_density, 2, raw = TRUE)          2 0.069360        21   0.065892
## aridity                                     1 0.038473        20   0.027420
## poly(shrub_density, 2, raw = TRUE):aridity  2 0.027420        18   0.000000
##                                             Pr(>Chi)    
## NULL                                                    
## poly(shrub_density, 2, raw = TRUE)         < 2.2e-16 ***
## aridity                                    < 2.2e-16 ***
## poly(shrub_density, 2, raw = TRUE):aridity < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m7 <- glm(animals ~  poly(shrub_density,2, raw = TRUE) * aridity, family = poisson, site_2023)
summary(m7)
## 
## Call:
## glm(formula = animals ~ poly(shrub_density, 2, raw = TRUE) * 
##     aridity, family = poisson, data = site_2023)
## 
## Coefficients:
##                                               Estimate Std. Error z value
## (Intercept)                                 417.057375  44.766454   9.316
## poly(shrub_density, 2, raw = TRUE)1         -34.363745   3.701733  -9.283
## poly(shrub_density, 2, raw = TRUE)2           0.703853   0.075745   9.292
## aridity                                     -36.376244   3.923030  -9.272
## poly(shrub_density, 2, raw = TRUE)1:aridity   3.022311   0.324374   9.317
## poly(shrub_density, 2, raw = TRUE)2:aridity  -0.061875   0.006638  -9.322
##                                             Pr(>|z|)    
## (Intercept)                                   <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)1           <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)2           <2e-16 ***
## aridity                                       <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)1:aridity   <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)2:aridity   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance:  2.7731e+02  on 21  degrees of freedom
## Residual deviance: -3.7303e-14  on 16  degrees of freedom
## AIC: 117.55
## 
## Number of Fisher Scoring iterations: 3
anova(m7, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: animals
## 
## Terms added sequentially (first to last)
## 
## 
##                                            Df Deviance Resid. Df Resid. Dev
## NULL                                                          21    277.313
## poly(shrub_density, 2, raw = TRUE)          2   65.371        19    211.942
## aridity                                     1  121.606        18     90.336
## poly(shrub_density, 2, raw = TRUE):aridity  2   90.336        16      0.000
##                                             Pr(>Chi)    
## NULL                                                    
## poly(shrub_density, 2, raw = TRUE)          6.38e-15 ***
## aridity                                    < 2.2e-16 ***
## poly(shrub_density, 2, raw = TRUE):aridity < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m8 <- glm(richness ~  poly(shrub_density,2, raw = TRUE) * aridity, family = poisson, site_2023)
summary(m8)
## 
## Call:
## glm(formula = richness ~ poly(shrub_density, 2, raw = TRUE) * 
##     aridity, family = poisson, data = site_2023)
## 
## Coefficients:
##                                              Estimate Std. Error z value
## (Intercept)                                 128.63472  111.86255   1.150
## poly(shrub_density, 2, raw = TRUE)1         -10.65865    9.24284  -1.153
## poly(shrub_density, 2, raw = TRUE)2           0.22049    0.18891   1.167
## aridity                                     -11.14405    9.80488  -1.137
## poly(shrub_density, 2, raw = TRUE)1:aridity   0.93610    0.81052   1.155
## poly(shrub_density, 2, raw = TRUE)2:aridity  -0.01935    0.01657  -1.168
##                                             Pr(>|z|)
## (Intercept)                                    0.250
## poly(shrub_density, 2, raw = TRUE)1            0.249
## poly(shrub_density, 2, raw = TRUE)2            0.243
## aridity                                        0.256
## poly(shrub_density, 2, raw = TRUE)1:aridity    0.248
## poly(shrub_density, 2, raw = TRUE)2:aridity    0.243
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 8.5698e+00  on 21  degrees of freedom
## Residual deviance: 8.8818e-16  on 16  degrees of freedom
## AIC: Inf
## 
## Number of Fisher Scoring iterations: 3
anova(m8, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: richness
## 
## Terms added sequentially (first to last)
## 
## 
##                                            Df Deviance Resid. Df Resid. Dev
## NULL                                                          21     8.5698
## poly(shrub_density, 2, raw = TRUE)          2   2.2783        19     6.2916
## aridity                                     1   3.8269        18     2.4646
## poly(shrub_density, 2, raw = TRUE):aridity  2   2.4646        16     0.0000
##                                            Pr(>Chi)  
## NULL                                                 
## poly(shrub_density, 2, raw = TRUE)          0.32009  
## aridity                                     0.05044 .
## poly(shrub_density, 2, raw = TRUE):aridity  0.29162  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m9 <- glm(evenness ~  poly(shrub_density,2, raw = TRUE) * aridity, family = gaussian, site_2023)
summary(m9)
## 
## Call:
## glm(formula = evenness ~ poly(shrub_density, 2, raw = TRUE) * 
##     aridity, family = gaussian, data = site_2023)
## 
## Coefficients:
##                                               Estimate Std. Error    t value
## (Intercept)                                  2.626e+01  4.945e-12  5.310e+12
## poly(shrub_density, 2, raw = TRUE)1         -2.167e+00  4.084e-13 -5.305e+12
## poly(shrub_density, 2, raw = TRUE)2          4.423e-02  8.343e-15  5.302e+12
## aridity                                     -2.288e+00  4.335e-13 -5.277e+12
## poly(shrub_density, 2, raw = TRUE)1:aridity  1.890e-01  3.583e-14  5.276e+12
## poly(shrub_density, 2, raw = TRUE)2:aridity -3.851e-03  7.319e-16 -5.262e+12
##                                             Pr(>|t|)    
## (Intercept)                                   <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)1           <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)2           <2e-16 ***
## aridity                                       <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)1:aridity   <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)2:aridity   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 2.626533e-28)
## 
##     Null deviance: 1.5420e-01  on 21  degrees of freedom
## Residual deviance: 4.2025e-27  on 16  degrees of freedom
## AIC: -1327.7
## 
## Number of Fisher Scoring iterations: 1
anova(m9, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: evenness
## 
## Terms added sequentially (first to last)
## 
## 
##                                            Df Deviance Resid. Df Resid. Dev
## NULL                                                          21   0.154197
## poly(shrub_density, 2, raw = TRUE)          2 0.137986        19   0.016211
## aridity                                     1 0.008285        18   0.007926
## poly(shrub_density, 2, raw = TRUE):aridity  2 0.007926        16   0.000000
##                                             Pr(>Chi)    
## NULL                                                    
## poly(shrub_density, 2, raw = TRUE)         < 2.2e-16 ***
## aridity                                    < 2.2e-16 ***
## poly(shrub_density, 2, raw = TRUE):aridity < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1